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ABSTRACT 


The problem of tracking submarine targets with passive 
sonobuoys 1S mathematically modelled. A comprehensive study 
is made of all the information available in the acoustic 
signals picked up by the sonobuoys and of the usefulness of 
this information in the estimation process. The presence of 
nonlinearities in the tracking model leads to the applica-= 
tion of nonlinear estimation theory. Bayes formulation con- 
cepts are applied to generate approximate solutions and 
filtering alaorithms, and the well known extended Kalman 
filter equations and higher order filtering algorithms are 


obtained from this approach. 


The concept of partitioning the measurements is 
presented and shown to brina advantages in computing effi- 
ciency and also, for nonlinear measurements, in tracking 
accuracy. A graphical interoretation of the action of Kalman 
filters is developed and provides insight into the %impore- 


tance of each variable in the filtering process. 


% 


Extensive simulations, desiqned to test the performance 
of the algorithms developed, are oresented in graphical form 


and analyzed. 
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SYMBOLS AND CONVENTIONS FOR TEXT 


l = Underlined small letters represent vectors,s thuS, 


x denotes the vector x. 


ec = Capital letters are used to represent matrices. 


3 = A caret indicates an estimated valuer e€.Ger x 
denotes an estimate of x. ‘ 
4G - Time points) are represented by t Ore when 


k 


between parentheses, simply by kr @eGer x(k) denotes the 


value of x at time te 


S = Estimates specify the time of estimation and the 
amount of information used, e.Ger x(kij) denotes the esti- 


mate of x at time t taking into account all the measurements 


k 
uo to and including et When k = 3} only one time point 1s 
indicated, thus, x(k) aenotes the estimate of x at time t 
taking into account all measurements up to and including the 


6 =2* Probability density functions arene pre Sern tie damety, 
the small letter pr 1-6€. Ply % 1 Y? denotes the conditional 
density of x given ys. When no confusion is oossibdle, this 15 


simplified to o(xiy). 


7 = The expectation operator 18S represented by Ef )3 


the variance by Var[{]}. 





a(k) 
a Ck) 
a p(k) 
a oC k } 
b(k) 


b (kK) 


Bk) 


expected value of the estimated state vector at Lo 
random rate of chanqe of target speed. 

random rate of change of target heading. 

random rate of chanae of frequency emitted. 

heading of target at time the 

bearing measurement by buoy i. 


second moment of the state estimate around expected 


value a(k). 


ee (tk) 

L 

kj 
e(k) 
e(k) 
f (k) 
ee k ) 
al 
Gk) 


HCk) 


Ek) 


average speed of sound. 


distance of target to buoy i. 


Kronecker deltas = 1 af k = jr = O 1f kK = je 


estimation error vector. 


expected value of the estimation error. 


rest frequency emitted by the target. 


frequency measured by buoy ji. 


gain matrix. 


observation matrix. 


comoression/expansion factor 


eS Ginait 1p error covariance matrix. 





QCk) 
R(k) 
s(k) 
Sk) 
0 Ck) 
o Ck) 


oCk) 


ie 
ij 


by buoys 


state excitation covariance matrix. 


measurement noise Cavariance matrix. 


speed of target. 


relative 


Standard 


standard 


standard 


time delay. 


time delay difference between the 


1 and j. 


time spent 


time soent 


deviation ofa 


velocity of target towards buoy 


deviation ofa. CL 


(k). 


b 


deviation ofa 


(k). 


f 


In prediction. 


in estimation. 


signals 


1. 


vector of q random measurement noise signals. 


vector of m random forcina 


inputs. 


state vector of dimension n. 


x*=position of target. 


Xx"=PDOSTtTION of buoy 


1. 


estimated state vector. 


[n€n+3)/2) -<dimensiona] 


10 


vector containing 


the 


received 


state 





variables and the distinct components of the matrix P,. 


y(k) y=position of target. 

yk) y*position of buoy i. 

zk) vector of q measurements. 

24 set of all measurements up to and including time ce 


Lae (2(0), z(1), oeeof 2(k)). 
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I. INTRODUCTION 


This work was oriented to the problem of optimally 
estimating characteristics and/or parameters of a certain 
class of nonlinears dynamics stochastic systems from 
observed output sequences which are noise corrupted none 


linear functions of the system states. 


In particular, attention was directed toc the problem of 
real-time, shorterange optimal localization and tracking o f 
a submarine target by passive acoustical means. The most 
accurate and reliable estimates of the target's parameters 
(position, headings speed) are sought by optimally process- 
ing the measurements obtained from the acoustic signals 
transmitted by the target itself and received by special 


sonobuoys. 


In the unclassified literature there are presently 
available methods of passive tracking of submarine targets. 
Reference [1] utilizes dopoler-shifted frequency measure- 
ments obtained from a group of sensors; (2) uses bearing and 
doppler-shifted frequency measurements ootained from one or 
more sonobduoys. These methods provide very good first 
approximations to the solution of the problem and many of 


their conceots and notation are used in this work. 


WS 





As a first step, an attemot was made to produce a 
comprehensive study of all the information available in the 
acoustic signals picked up by the sonobuoys and its useful- 


ness in the estimation porocess. 


The estimation algorithms were developed taking into 
account the advantages of processing Information as soon as 
it becomes available. Flexibility to move, remove’ and 
include sensors and measurements at any time was obtained. 
Constraints were included to account for limitations) of 
practical, inexpensive and presently available sonobuoys and 


computing equipment. 


In Chapter II the problem is described in mathematical 
terms using state space techniques and the initial assump- 
tions and constraints are given. The model used for the 
target 1s similar to the model presented in [2]+- but with a 
different set of states. The measurements obtained from the 
signals received by the sonobuoys are characterized and 
their functional relationships to the parameters of the taro 


get are analyzed. 


Chapter III discusses the general theoretical solution 
for the problem and the difficulties of implementing it. 
Approximate solutions are then sought and processing equa- 
tions are developed. The special vectors, matrices and 
Relations characteristic of the tracking problem are 


prepared for the application of the filtering equations. 


14 





In Chapter IV the concept of partitioned measurements 
1S Introduced and analyzed. The advantages in accuracy, come 


puting efficiency and filter flexibility are stressed. 


Practical and graphical interpretations of the estima- 
tion process help in visualizing nonlinear problems and 
motivate the adoption of iterative techniques. Increased 
accuracy and possible converaence improvements are shown to 


result from this aporoach. 


Chapter V introduces the interactive computer program 
that allowed the apolication of the ideas describea in the 
previous chapters to the tracking problem. Selected simula- 


tions are discussed and presented in graphical form. 


This research covers the subject in @ much more 
comprehensive way than earlier works, [1], [2], [3]. The 
results obtained by using the methods and ideas collected 
and developed in this study show many ways of improvina the 
tracking accuracy, of increasing comouting efficiency, of 
reducing divergence problems and of obtaining faster convere- 


gence. 


The importance of having the buoys placed in adequate 
positions, of having an adequate model for the target, and 
of having freauentr varied and precise measurements is 


Stressed. 


15 





The final chapter Summarizes the results of this invese= 
tigation and presents the conclustons and suggestions for 


further study. 
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II. PROBLEM FORMULATION 


A. THE ESTIMATION PROBLEM 


A general nonlinear, stochastics dynamic system observed 
by a set of nonlinearsr nonstationary measurements can be 
represented, for our ourposeser in state=space discrete form 


by the equations 


xCk+1) = FOxCkI owlkdet, s kee] 


z(k) = A(x(k)rv(k) rt, ) (2.2) 


where x(k) is the state vector of dimension n, 
wk) 1s the vector of the m random forcing inputs, 
z(k) 18 the measurement vector of dimension a, 
v(k) is the vector of p random measurement noises, 


f( ) and h(€ ) are general vector functions. 


The general estimation oroblem consists of 
“-= knowing the statistical characteristics of wk) and 
v(k) and having a statistical representation of the initial 
conditions x(0) of the system, 
--- processing the observations 2(0), z(1l), «seoor 26j) 
in a statistically optimal way to obtain the best estimate, 


iM some senses of the system state x(k) at time the 
If k > j the processing is called predictions, if k < j 


it 1s called smoothina and if k = j it is called filtering. 


a 





In this work some departures from the completely general 
nonlinear case were taken. The first assumption is that the 
measurement noise 1s additive and that the system equations 
are linear with respect to the random input vector we Equa- 


tions (2.1) and (2.2) can then be written as 


xCk+ld f(xkirt ety + 9 OSCE EA EAN, Ae? 


Cano 
z(k+t1) = h(x€k)yt, ) + v(k) (2.4) 


mere f€ J) and ht ) represent general vector functions and 


gf ) a general matrix of functions. 


Ssecondliyr the measurement noise and random forcing input 
are assumed to be uncorrelated in time, zero mean, discrete 
Gaussian sequences, independent of one another and Indepen- 


dent of the initial condition of the state of the system, 


ivel es 
Elw(k)] = 0 ey Ce) =a C2 as) 
lee . Mee. re 
Elwlk)woCj)d)] = Qk) 6, Elvtk)veli)] = R Games . 
— kj ah, kj 
apa . T = - AL z 
Elwlk)v (j)) = 0 Efx(O0Jw (k)}] = 0 Efv(k)}x (O))] = 0 


| 
7 


-" 
Coe = 
, 
J Oiftk #j 





Be IME PASSIVE TRACKING PROGIEEM 


Many practical situations can be formulated as estima-~ 
tion problems and characterized by the equations described 
above. The Situation that motivated this study assumes. a 
submarine target following (most of the time) an approxic 
mately constant speed and constant heading path in ae field 


of passive sonobuoys. 


The target unintentionally emits acoustic signals which 
are picked up by hydrophones. The transduced signals are 
sent up by wire to the buoys and then frequency=modulate VHF 
carriers which are transmitted by the buoys and received at 
anearby ship or aircraft. After the recovery of the sianal, 
data processing equipment generates measurement values which 


contain numerical information about the target parameters. 


The first step in the analysis of the signals collected 
by the sonobduoys is the attempt to detect the presence of a 


real target. 


After a target is detected, frequency spectrum informae- 
tion is added to intelligence data in an attempt to classify 
it. Other measurements and information from any other source 
are also used to obtain aefirst approximation to the 
target's parameters (position, heading, speed), generally 


through the use of least mean square techniques [1],[3]. 


The third phase of the processSer if a wartime condition 


doesn't exist, is the trackina of the target. That iss the 


NS 





determination of its eath while reducing, if possible, the 
initial uncertainties about its parameters, in real-time, 
through the optimal processing of the sparse information 
provided by the sonobuoys. The passive characteristic causes 
severe limitations but is necessary to avoid revealing to 


the target the fact that it 1s being observed. 


This investigation was devoted to the tracking phase and 
it was assumed that the main characteristics of the buoys as 


well as their positions are known. 


C. THE MODEL OF THE TARGET 


The problem isvr of courses three=dimensional. Neverthe= 
less the measurements are not a direct function of the 
target's denoth but of the difference in depth between the 
target and the hydrophones. Anticipating the accuracy of the 
measurements and the precision of the data processing equip-= 
ment and algorithms, it is justifiable to consider only two 
dimensions initiallys, adding depth as the third dimension 


whenever necessary and with no conceptual aifficulty. 


The target's path is frequently one with nearly constant 
soeed and heading, disturbed by currents and random thrust 
and control variations. Intentional maneuvers which have no 


evasive purposes are normally simole and smooth. 


A basic olant havino as states two position coordinates 
(x anc yoy heading (bo) ond speed (s) of the target, descri a= 


INQ a Constant velocity paths, as shown in Fig 1, seems 


20 





adequate. 


Random forcing inputs should be considered to account 
for noise processes and imperfections of the model, such as 


target maneuvers. 


As suggested by (2) these random forcing functions can 
be approximately represented by two independent zerormean, 
piecewise=constant random rates of changer a | and Oe act= 


ing on the speed and heading of the target. 


The formulation of this basic plant in discrete  state= 
space form, similar to Equation (c.3)+r can ve obtained in 


the following way 


t 


k+1 ; 
= = = = - k 
s(k+1) s(k) a(t) dt WG ae th? a ¢ ) 
“k 
“etd : 
b(k+1) - bk) = = : ° k) 
) Ck) : a(t) dt ee i af 
k 


ti 
it? 


Sell 
s(t) cos S(t jet 
t 


pe . 


x(k+ti) = x(k) 


{s(k) + (t - ti Ja (kd) cos{b(k) + 


Gk 


+ (t ° tp) Ck] dt = 


Geaal . 
s(k) cos b(k) dt + 
t 


Sst 
+ (t - ti Ja. Ck) cos b(k) dt = 


Bie 


ce 
[~“ - t da, Ck) sk) sin b(k) dt + ws. = 


e 


k 
(tia - eu s(k) cos bk) + 


a) 





target position 
ayie teal typ 





Figure 1 : Geometry of target and sensors. 


sage 





1 Z 
—— - k k - 
5 Oy t FCS eos b Uk) 
- oe wlse s(k) sin b(k))] + HOT (2560 
where HOT represents Higher-Order Terms. 
In the same way, 
Cet 
Vat i) = y(k) = s(t) sin b(t) dt = 
. t 
k 
= ad + 
Coe t s(k) sin bC(k) 
1 2 ; . 
+ — t - kk bk + 
5 ( al Bo bo Visine Ge) 
+ a CK) s(k) cos b(k)) + HOT Cen) 
I f (tray - t. Jak) and Cty - Cor (kK) are suffi- 
ciently small so that the higher order terms in Equations 
(2.6) and (2.7) can be neglected, one hase in vector form 
x(kt+1) = f(xCk)et oti) + SMES OIG Da) oll se 
where 
x(k) 
x(k) = y(k) (2.8) 
s(k) 
b(k) 
a (Ck) 
Ss 
w(k) = (2.9) 
kk 
as | 
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f k t)) = Gece) 
See See Be 
yl ee cae - t ) s(€k) cos bk) 
k+1 k 
y(k) + (Ct - t ) s€k) sin bk) 
k+1 k 
s(k) 
b(k) 
and 
Cle jig t t ) = C2n1)) 
Ss rans 
il Z 1 2 
a(t yay - t.) cos b(k) 5 (tay t s(k) sin otk) 
It ae i Z 
7 S| t) sin oO(k) Geral ee ee cos b(k) 
eo aa = t,) Q 
0 cen - t) 
When frequency measurements are to be orocessed 


directly, as shown in Section I[1I,D,1,0, an additional state 
variable must be included "= the rest frequency oe emitted 


Dy the target. 


This frequency is assumed to be approximately constant 
with ae$smalj) random disturbance Q ¢ that 1s assumed to be 


zero-mean, piecewise"constant, and independent ofa. and Oy - 


The expanded state variable formulation iS aiven Dy 
x(k) 


x§(k) = 2--- (2.12) 


f(k) 
0 


—/ 
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wk) 


w°(k) = ---- (2.13) 


a Ck ) 


me Okt!) = fF xe Ck) ety, rt) + 


if 


a g “(x° (k)-pt 


Wee ee Cama 


where 
flxCkir ti rt) 
e = nn aR erica, Sk 
at Cx Ck)rti i, ety) = (Care tos 
f (k) 
O 
| 
e,_e = 
g Cx CkJet et) ee to (An he 
| 
0 0 ! - 
Nee eae 


c 


(2.16) 


D. MEASUREMENTS 


The songbuoys can perform one or both of the following 


tasks: 
“"-= pick up underwater sound signals 
@“-- indicate the approximate direction of the vectorial! 


Sum of the received acoustic sianals. 


a5 





The target emits a siaqnal r(t) which 1s distorted = and 
modulated by the propagation medium between the source and 
the hydrophone, and by extraneous sound sources present mn 


the ocean. 


When all of these disturbances are of relatively mild 
Strenaqths the signal received by a stationary buoy i will be 
approximately an attenuated and delayed version of the pre-e 


viously emitted siqnals ieee, 


Hie 


r (t) a (t).rCtme 1. (t)) (2.17) 
1 aL 1 


where a) is an attenuation factor and tT, Ct) is the time 


delay given by 
Tact) ==—deXtt —= 7 (i) Sic C2.185 
ah 51 i 


that is, the distance d_ between the target and the buoy at 
i 
the time of emission divided oy the average velocity of 


sound in the medium. 


Since the target velocity 1s much smaller than the sound 


velocity, for small distances and delays one has 


~ 


Cad 


ect tT Act. = Ot) + s (t)er (t) (2.19) 
1°: i aL i a1 


where s (t) is the relative speed of the target toward buoy 
i 


1 aS Shown in Fig ee. 
Thuss from Equations (2.18) and (2.19) one has 


t.(t) = da.(t) J Ce = § (t)) (2.20) 
a aE 1 
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SS position of target 
at time t. 
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d(t~-t(t)) 
i al 





Figure 2 : Target-sensors geometry - distances 
and relative velocities. 


bad 





and the received siaqnal 18 given by 


eo (t) = a, C(t) elt. = dete) Abem= s. Geo Game) ) 
aL z aS aU 1 


If this signal iS recorded from t = e to t = ¢t 


one has recorded 
ttt | i seca Otome tt" ermrter t vt! = 
io O 
=Go2Ct. tot i/ Cee = Sots. Steamy 
al O 1 O 


ass t' < T 


If during this oeriod T the relative target 


towards buoy i doesn't chanaqe significantly, 
Sat weet) s= Ss Ct. 3 
. | re) 


secure) ="dmet ) - t'.s (t ) 
A ls O 1 O 


speed 


re (te + t') Fa (Ct + t').r(t = dace .= ss Cu mat 
Lao O 1 Oo j 


1 


tti.ft t+ s(t /(e - s Ce 2)]] 
at O al O 


me Gt, tint!” i= on (te ate Cena (Genin Gate ))/ (cules CUNO amt 
O 1 


LE iO aE O O L 


t Cate 7tc = "Ss, Cie) Canen) 
ali ee. 


In Equation (2e.ec), the second term inside the brackets 


1S a fixed time delay. The third term is a variable time 


Gelay and can be seen as aA comoression/expansion term, 
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ing the different value of the variable time for the origi=- 


nal and the received signals. 


If the emitted signal had a single frequency component Wo 


one would have 


p(t Sr t' ) Ss Cte tC COS tHE Gte tt ID 
ibe Ke: igmo io 


where 


tt 
the 


cos(w (t + t')) 
1 


cos(¢ + w et’) 
oO 1 


eo) 1 


Clos) tw! foreman We Ge tGcem i (Culm aS.) ( ba) at 
O O O L O 1 O 


Poesia) Goa CaS. (te) 7) = 
O L O 


Cos Gow tow .c.ctmy fe = s (t ))) 
10 O io 


and the received frequency would he 


<z 
1 


Wie G i, UC CamorS 
a 0 1 


eee tf oC US CG Se.) (2.23) 
The siqnals collected by varicus sonobuoys are sent to a 


ship or aircraft normally equipped with equipment and skills 


Bowextract the information described below. 
1s Isolated Buoys 


By processing the signals collected by each buoy 


aloner one can obtain the following noisy measurements: 
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ae Bearing Indication 


Bearing 18S aQgiven directly by the buoys) or 


preprocessed (filtered) to reduce noise effects. 


Since the bearing is obtained from ae vectorial 
Sums it can be reasonably accurate only when the acoustic 
noise 1S approximately omnidirectional or much weaker than 
the target siqnal. For some sensors this vectorial sum can 
be limited to a selectable frequency band, thereby easina 


the noise problem. 


The relationships between a bearing indication 
obtained by a buoy 1 and the states of the target 1S shown 
in Fig 3 and given by 


Ck) ) + vk) 


b.(k) = arctan[{y(k) = y.(k)I/Ex(k) = x 
1 1 1 


(2.24) 


where v(k) must account for all the noise aisturbing this 
measurement, including acoustic noiser transducer inaccura= 
cies, preprocessing noise and errors due to the finite speed 
of the sound in the water. This last factor 1S caused by the 
time delay between the emission of a sound by the target and 
reception at the hydrophone. At the end of this period the 
target and the hydrophone have slightly different relative 


positions. 


t 


According to (2] tyoical accuracies of + 5 
degrees are common for strong sianals and inexpensive DIFAR 


bDuOYSr and as cited by [3] tnaccuracies down to +t l degree 
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Bh 


Figure 3 


x (k) x (k) 


b{k) = Dts 
me 


: Bearing measurement. 


ail 








or less can he obtained with arrays of hydrophones and 


proper preprocessing. 
aye Frequency Indication 


By analyzing records of lenath T seconds of the 
acoustic signals picked up by the sonobuoys with Fast 
Fourier Transform algorithms, one can obtain approximate 


frequency spectrum descriotions with a resolution of 1/T 


Hert Ze 


By detectings measuring and possibly tracking 
{4] some of the strongest frequencies contained in the sig- 
nal», very useful data about the target parameters 1s 
obtained since the received frequencies are functions of the 
originally emitted frequencies, the soneed,s, the heading and 


the position of the taraet, as shown in the relationshio 


f (k) fai) ec GG sok) ) pat OAC ka Gace) 
di O iL 


where 


Ss Ck) “s(k).cosib(k) = b. (k)) C2..c60) 
1 


is the relative velocity toward buoy ir as~shown in Fig 3. 


The term v(k) now has to account also for the 
errors introduced in the various paths of the signal from 
the target to the ship or aircrafts and for variations in 


target~sensor geometries during the record time I. 


3e 





As suqgested in (1), (2) and (3), inaccuracies 
ranging between 0.5 and 0.01 Hertz can be obtained in most 
situations for values of un of the order of hundreds of 


Hertz. 
Ce Signal Strength Indication 


The strength of the signals varies not only with 
the strength of the originally emitted sound but also with 
the distance between the target and the buoy. The signal 
strength also depends on the random aspects and the direc 
tivities of the target, the hydrophones,y, and of the 
transmitting and receiving VHF antennas. Random absorotion 
and scattering in the various media from the target to the 
processing equipments may also have a significant influence 


on the signal presented to the processors. 


A mathematical formulation of all these effects 
would be extremely difficult task and would obscure the 
desired distance information. For this reason the informa- 
tion contained in the strength of the collected signals, 
which can be used by an experienced and skillea audio opera-= 


torr, was not incorporated in this study. 


Attenuation factors are thrown away by setting 
the strength of the received signals at a good working 


level. 
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2% Two or More Buoys 


Considering the signals received by two or more 
buoys, indications can be obtained which may be used in 
addition tors or instead ofr, the measurements available from 


a single buoy. 
ae Frequency Difference Indication 


To process frequency measurements obtained from 
one buoy requires inclusion of the rest frequency as one of 
the states. This is the case because this measurement is 
very sensitive to eee in the assumed rest frequency, 


TeCe, 


Af. = Mf. Cf SCC = Ss.) = Af 
a: O a O 


Having the frequencies received by two buoys and 


taking their difference one gets 


Wee wtio = fevce7 Ce = SS.) = £..6¢7 Fe = Ss.) = 
J O 1 O J 


fe. US. = —S.9 07 Tle = Ss.) (eas Sila] 
O 1 j Be i] 


We 


O 


{eS asi . (2367) 
4 4 


The sensitivity of this difference to errors in 


the assumption of the rest frequency is given by 


Cc 
oo 
hn 
4. 
) 
—h 
Se? 
vt 


Aft=.(s .= s.)/c 
O al j 
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or is reduced by avery large factor since c is of the order 


o f 1500 m/s and Cs. - S27 is normally very small. 


This way, if one has a reasonable approximation 
for fo , instead of processing two measurements ina five= 
state plant one can process only one frequency-difference 
measurement ina four=state olant with greatly reduced com= 


puting time and hopefully not a detectadle reduction in 


accuracy. 


The variance of the errors in each measurement 
must be added to give the variance of the error in the fre= 


quency difference measurement. 
be. Frequency Ratio Indication 


If one divides the frequencies measured at two 
sonoouoySs,r a new relationship 18 obtained that retains the 
information on positions, heading and speed of the target, 


but. 1S indeoendent of the rest frequency emitted 


fee?’ tC elt nc co = Ss) la ilcen is) 0/7 tt ..c)) = 
J O i J fe) 


=e C 4S) WNC ce= ‘Ss 9 CO eo 
j il 
or 
efoto =.) CS ee er) Cece ele See Cart (s/c)? +) Seren = 
1 J 1 J a i 
= + (sees 07 -c (2220) 
1 J 
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As with the frequency difference indication, 
this measurement can be used in a four=state plant instead 


of using two frequency measurements ina five=state plant. 


If this ratio is generated by the division of 
the frequencies obtained from the analysis of each of the 
Signals, the characteristics of the frequencysratio measure= 
ment noise is complicated and state dependent. If however 
this relation 1s obtained as described in the following dis- 


cussion, this uncomfortable situation can be avoided. 
Ce Time Delay Indication 


As shown at the beginning of this section on 
measurements, the signal collected by a sonobuoy 1S a noisy, 
delayed, attenuated and compressed/expanded reproduction of 
the original signal emitted by the target. The delay is due 
to the finite time the signal takes to reach the buoy” and 
the compression/expansion is caused by the variation of this 


delay with time (doppler effect). 


The signals received by two sonobuoys have adife= 
ferent noise contributions and different delay and 
compression/expansion factors even after their strengths are 


equalized. 


Supoose one recorded the signals received oby 
buoys 1 and } irteigh 8 = t. to t = t + T. The signals are 


related to the original signal as shown in Equation (2e.ce), 


Te@er 
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r.(t. 4p!) “Samer + ter lt) = Sancta Ce 
i= 56 i 6 O io io 


tis¢wt 7 Ge =-.s2) (ta) } (2.350) 
Al O 


p(t + Ce. Che) (tac WG) 7 ce = ct ee 
j- o jo fe) lee j ae 


tac en 7 (6 = >s a(t.) ) (2.351) 
Jj O 


If one now amplifies these two signals to a good 


working level, correlate them by evaluating 


T 
pC = | ri(t + t')er C(t + t' +7 )edt! (ase) 
i eo j Oo 


O 


and look for the value of t that maximizes o (1) one’ finds 


thats: 


C1) -1f *s, (t>) = s(t ) ands, consequently, the fre- 
ao q O 
quency shifts are about the same in both signals, then r (Ct 
io 
+ t') will be simoly an approximate delayed or advanced 
replica of r.(t + t'),» as can be seen from the above equa- 
J O 
hilons . 


In this case the maximum of o(t) is easy to 


find and corresponds to the value ae at which 


t = d (te y/le = 5 (t )) + calt* + 4 J/€e = gs C(t )) = 
fe) i 6 a Oo 3h {o! 


max 


= fe Ot tema S Ch) whecet 7 CCe = pSemetn ) 
O j 0 J oo ja 8 


Solving for ft yielas 
max 


T = (dct) eager )) 4 € Ceeoo9 
max io iO 


or 





(2) if s.€t_) FJus.(t@)> pp (qe) deecmenctaehave wa 

Cea, je 
simply=determined maximum value because the two signals have 
frequency components whose phases are shifted by different 


amounts. 


If ones now change the time scale of one of 


the signals, say pus + t’), one has 


r'(t') = r,t + met’) = 
le J O 


Um: 


a, or ({t =“d (tay) / (c = ase(t 0) «Ct 
J O J O J O 


tcemet*/(c = s.(t_))] (2.34) 
J O 


Comparing this equation with Equation (2.30) 


one sees that if m can be found such that 
meer (Cem 6 (t_)) = c/Cce = s (t )) 
jo io 


or 


me= (Cc = Suwlt®)] / fe =. s(t 3) (2555) 
jo io 


one again has case (1) and the correlation function wil) 


give a maximum at aporoximately 


~~! 


t = {d.(t ) © d.(t )) /e (2.36) 
max 1 oO J O 

It is interesting to note that Eauation (2.35) 

is the same as Eauation (2e.2c8). The compression/expansion 
factor is thus the ratio between the frequencies received by 


buoys 1 and je. 
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The p (tem) functions if the two recorded sig- 
nals were broadband, of long duration and with high S/N 


ratios, would appear as shown in Figo 4. 


Low signal to noise ratios could hide the rea) 
maximum and create false ones; narrowband or coherent sige 
nals would present ambiguities with the creation of many 


close maxima. 


By considering the availability of equipments 
and techniques to obtain the proper elimination of the 
difference in frequency shifts 1n the recorded signals , 
this study assumed the possibility of having time delay 
measurements available to provide information about the tare 


get position accordina to the relationshio 


Peck) = td tk) md Ck)ae 7 ¢ + vik) Ud aaa) 
1j 1 a 


where v(k) must account for the errors caused by the dif-= 
ferent noise contributions tin each signal, for the accuracy 
and orecision of the cross=correlation algorithm or device, 
anacd for variations in target-sensors aeometry during the 


recording of the signals. 


Inaccuracies anywhere from 0.5 sec down to a few 


milliseconds seem very reasonable to exoect in many cases. 
d. Signal Strength Difference 


For the same reasons already explained in 


My Drlsc the information contained in the difference of 
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(time delay) 


(compression factor) 


Figure 4 : Correlation function. 
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strength of the collected signals, while of great use for an 
experienced and skilled audio operator, was not incorporated 


in this study. 
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TIT. THEORETICAL SOLUTION 


A, RECURSIVE BAYES FORMULATION 


The maximum amount of information at time tee in a pro- 
babilistic sense, is concentrated in the conditional density 
function of the states given all the observations up to this 


time, (x(k) $2"), where as Represents the set i2(k) 7 ez Cee 


eee f z(0)}. 


By specifying a cost function one determines how this 
information is used to obtain the best estimate [5]. For 
example, the minimum variance estimate is the mean of the 


eonaitiona! density. 


From the assumptions made in Section II,A, the random 
processes x(k) and {x(€k),-z(k)} are Markov and the condi- 
tional density can be written in recursive form using Bayes’ 


Lawe as is shown in (5), (6), (7), (8). 


p(x(k) $25) = ek) ox (k) LZ wp (2 (kd Px (kd) (Ga 
where 

pCx(k EZ) = | px Cke1)FZ4 p(x (kd bx (ko) ) dx (ke) 

: CSesap 
and 
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17 eC ee ee 


= | P(x Ck iZ S| ote Ge) x(k) Ie dxk) (3.1b) 


Given the initial condition pier Oza) = p(x(0)) and the 
probability law of the system random input and measurement 
noises, plwlk)) and p(v(k)J,one can proceed to evaluate these 
equations and the friterinag problem can be regarded as hav-~ 


ing been solved when peck) iz~) can be determined for all k. 


For any Situation different from the one where the sys- 

tem and measurement equations are linear and the apriori 
: : ; ; k : 

distributions are Gaussian, the density p(x(k)iZ° ) is not 


Gaussian and generally cannot be determined in closed form. 


In the problem studied here, the conditional densities 
in Equation (3.1) are nonlinear functions of the states and 
the measurements. Linearizationse, Taylor expansions) and 
Gaussian approximations will be used, however, aS an 
engineering comoromise, 1n order to avoid complex numerical 
integrations in nm dimensionse ThisS approach was chosen tak-~ 
Ng INtO account the limitations on processing equioment 
available for the practical solution of the addressed track= 


ing problem. 
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B. ONE-STEP PREDICTION 


Let us assume that at time . , Guring a filtering proo 
cesS, an approximately Gaussian density o(x(k-1)!1Z*) has 
been obtained with mean value up (k-1) and covariance matrix 


Bik = | 


If one chooses the mean of this conditional density as 


the estimate at time eae 
e(k-1) = Etx(k-1)iZ24 = (ket) 


the estimation error at time k=l has the same probability 


density as o(x(k-1) $2574 but with zero meaner 1.@.e,r 
Ele(k-1)} = EU(x(ke1) = x(k-1))42K74 = 0 (3.2) 
Ele(k=1tet(k=1)) = P(k-1) = M(k-1) (505) 


At time t anew set of measurements arrives and a pred 
iction~ at that time must first be obtained based on the 


measurements u to ¢t - 
E eon 


mibie dynamics of the system are described by Equation 
moe and it can be seen thats even if p(x(ke1) 4; Ze) were 
really Gaussian, the density soa) Ae would not normally 
be because of the nonlinearities involved. Nevertheless it 
will be assumed that a reasonably close approximation to a 
Gaussian density exists and proceed to find aporoximate 
values for the first two moments of the conditional one-step 


pre=-diction density Stn Oe) ieee se 
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The first step is towexpanc Equation (2.5) 











iv) @ way lor 
series about the estimate EBS) « The dependence on t. and t 1 
is dropped for simplicity. Note that z*-1;. known and was 
used jin determining x(kel). 

x(k) = f£0xke1)) + glx kel))wlke-1) = 

= f(x(ke1)) + us o(x(k-1) = x(k-1)) + 2... + 
='K(kK-1) 
+ a(x(k-1)).w(k-1) + ... = 
n 
= Ce es af 7 
= f(x(k-1)) a a re he Lb || 
x(k-1) 
n on F 
+ ) ) ee ee (kel).e. (kl) + «es. of 
joe OX, OX, i j 
: &(k-1) 
Pe GCx Ck iwc 1) i 
5 122 
- = ew (kei) | .e(k-1) + ... (54) 
ge Ns | ol a af - 
pei) 


> ae 
where i 


and e, are individual components of x and (x - x); 
Q. 18 a column vector of g + and w. 1S a Component of we 
i = ll 


1. Linear Approximation 


In this case only the terms in the 


expansion which 
are linear in Gest 1h) and w(k-1) are considered, 1.@e, 
x(k) = f(x(k91)) = o(kel).e(ket) + 
t g(x(k-1)) ewlk=1) (3.5) 
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where 


o(kel) = f(x(k-1)rt, rt ) (3;56)) 


Kea 


In |~ 


From Equation (3.5) the mean value of p(x (kar oe 


which 18 chosen as our predicted values is 


Repke) = EUxCkKIIZ A = #eKRCKR1)) - 


- (kel) .Ele(k-1)) + g(i(k=1)).Etw(k-1)) 


pipe xk=1) and mw(k=1i) have zero means, therefore, 


Rakim 1) SOS IO 1g? COR 


Thus the one=step apeqiebias error has the mean 
Ele(kik-1)} = E((k(kikel) = x(k} = 
= Elglk=1).e(k=1) = g(x kel) ).wlkel)) = 0 
The prediction error covariance matrix, which is 


also the covariance matrix of Cade: ale iS, using simpli 


fied notation, / 


P(kike1) = Ele(ktkeLe (kik-1)) = 


e [tse ete - a (Fe 1) ou C2191 Go ee)7] = 


Elgeece 6 > Gres eg) = Gewee Bg™ + gewewT gl) 
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From the assumptions) about the independence’ and 


discrete white noise characteristics of wr the terms 


g (ke1) -Ele(ke1) wi (k-1)) a (R(kH-1)) = 0 
g(x(k-1)) Elw(ke-1). e (kel)].g¢ (kei) = 0 
and 
P(kike1) = > (k-1).P(k-1). (kat) ‘ 
t g(x(ko1)).QCk=1).a7 (x(k=1)) (38) 


2. First Order Aoproximation 


In this develooment one retains al] terms in Equa- 


tion (3.4) up to first order partial derivatives, 1.@e, 
x(k) = f€x(k-1)) = ofkel).e(k-1) + 


+ alx(k=1))ew(kel) - 
mM 
- J A, Chl) ew, Ck91) 6 e(k=t) (3.9) 
. a L io 
where 


=e tae 
Ap (kr1) 7 52 ay (aCe et et) (ic 0) 


As in the last sections the predicted value is 
x(kiket) = f(xCkel)) =9 (kel) .Ele(k-1)) + 


+ g(x(k-1)).ECwlk-1)) - 
at 


A .Elw. (k-1).e(k-1)) 
l 1 1 ae 


i 39 
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From the oassumptVons oebout witk= |) sncmee Ck =) 7 sean 


prediction 1s 
SSS SCS CRS 
the same as In Equation (3.7). 


The one-step prediction error has zero mean and 


covariance matrix given by, in simplified notation, 


P(k{k=-1) = E(t 1) een - g(x(kol)).wlke1) + 


be 
+ ) beet on, eed) eC) Get 
i i ~ 
i=l 
a TA 
= Elg-e.e! .of ee eeewi.gl + ) gele€ oA ew = 
_ oe an a ae 
m 
- Guess. 5) + gewewt.gl - J gewee .Al iw + 
ae ina = =e alo i at 
m m 
ae ) A sece ag ow, - ) A se ew ogi ew + 
i=l = - it—JP i:-7 a 4 
mn m 
+ y } ne. he ow J 
a ae L:_=->-- J 
1=1 j=1 


From the assumptions about wer one has that 


Ele(k-1).w (k-1)1 = 0 


Ele(k=1).e° (kl) ow, (k-1)) = Efe(k-1).e' (k-1)). 


Elw (k-1)) = 0 
ale 


Elw(kel)ew (k-1)). 


Etwlk-1).el(kol)ew. (k-1)) 
; 1 1 


Ete'(k=!)] = 0 
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EleCke1) se (k-1) em, (ko 2) ow, (k-1)] = 


Ele(kr1) sec (kr1)) E(w, (kot) ew, (k-09] = 


P(k=-1) e q.. (kel) 
alee 


where oT eae. 1S the jrlth term of Q(kel). 
Thus, 


P(kik-1) = (k=1).P(k=1) 6 (kel) + 


tig tk-1)).0(k=1}.q° (x(k=1)) + 


m m = 
+ } } A .P(k91).A. wa., (kO1) (3) 
=1 j=l J Jil 
If the @ matrix is diagonal, Sle = 0 for j # 1 and 


P(kikel) = o(ke1).P(k-1).o (kel) + 


+ g(x(k-1)).Q(k91) .a! (X(k=1)) + 
m 


+ JY A, PC ket) .AT Ca. (k=1) (3.12) 
4= iL 1 a LL 


Sve Other Anoproximations 


Depending on how many terms of Equation (3.4) are 
considered, other different prediction values and prediction 
error covariance matrices are obtained. Higher moments of 
the conditional densities would be required with increased 


complexity and computational burden. 
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Gye UP=-DATING THE ESTIMATE 


Let us assume that we have an approximately Gaussian 
density p(x(k) $784) whose mean value 18 the prediction 
x(k}k=1) and whose second central moment is the prediction 


error covariance matrix P(kjke-1l). 


A new set of measurements 18 givens Consisting of non= 


linear noisy observations of the states, represented by the 


relationship 


z(k) = hxCk) +t ) Feayvrek:) 


MES iS Mow mMecessary to process 2(k} and extract the 
information it brings. A new approximate Gaussian density 
p(x(k) $25) is assumed to be generated whose moments are the 
new estimate x(k) and the estimation error covariance matrix 


P(k). 
This new density is given by Equation (3.1) or by 


p(x(k)$2") = p(zCkedixe)) CeCe IZ 7 pCzce yp Zk-} 


aye tS) 


k- 
p(x(k) 12 5 1S assumed to be approximately Gaussian with 


moments x(kike-1) and P(kik-1). 


Expanding the measurement equation around x(kik=1) one 


has 


; h : 
ZC) = ACRE EK1)) + SB) ix Ce) © Rebel) # oe. # vk) 


x(k k=) 
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Retaining Only the linear terms, one has 


2(k) = AOx(kikel) # Heke i(k) = x(kEkKe1D) + ek) 


(3.14) 


where 


H(k) HCC k= It) (3.15) 


I | 


From Equation (3.14) one obtains 
Etz(k)iZ4 = hex(k}k-1)) 
and it is easily shown that 


Var (z(kK) $204 = HC kK) PC kEkeL) oH Ck) + Rk) 


Also from Equation (3.14), 


Bleak) ix (kj) = hidaCkiko1)) 4 mk). ix (k) = x Ck k= ton 
Var (2(k)ix(k)] = R(k) 


If one now assumes that M2 Canin, and p(zCk)ix(k)) are 
approximately Gaussian densities with the moments given 


above, then Equation (3.13) becomes 


oe 


P P 1/2 
p(x(k)$Z%) = [HC kK) .PCKE KEL) CHT (Kk) # ROK] / 
1/2 Wiz 
/ (C2nP/A(Peki ket) f .fRCk)] 1. exp [B) 


where 


a4 





Bete tone e fix i BODIE eet Cee] : 
' Ec - (x(k ket) = Hk) eix(k) = kCetketdd7, 
Tote | Z 


- [rece - h(RCkike1))I7. 


CS loD 


Simplifying the notation, this exponent can be reare 


ranged in the following way 
(x = 2) PL =) +e fz 2 he) - Hee - DIT RAL 


fz = h(s)i-. (H.PoHee RT iz - hO2)T = 


=x) Pa tx x) + Cx © xytHt etc, = %) - 


if 


z= nex) Rent x =) 2 =e Cee ln RE Se neeeniaees 


poe en Ce eRe me ee 


Seen Ce inen «4+ Rime eh cL Tee 


= 


Ge ae) ep = ee RiGee? Shee 


re es oe 
= {eo = Ce) Ren ee 


1 


th eee are oe ne ey 


Se 





Using the matrix inversion lemma [9], 


whe cHPHT ¢tey t= R HCP HT RAHStHT RS 


Since R and P are symmetric matrices, 
RHP + HR HSH R= tee ts HA 
where 


pte Hen el 


Be 
ot 


and the exponent of Eauation (3.16) becomes 
{x = oy Ges He rly Cy - x) - 
Sa Ole rab a ee 


OP i Ce Porn Bae 


~ (CP + He lnytH tet, = nok) cpmte ATH CL. 


amach 1S a Quadratic form that can be expressed as 


E -< - Pte HRI RT - cay] tp he Wet, a 


And finally the conditional density of Eaquation (3.16) 


can be put in the form 


N lk 


o(x(k)!Z") = al(k).exo P x(k) - Fook. C003) 


where, by definition, 
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x(k) = xCkiket) # PCKIHTCKILR HK). (2k) = NCR KEK) et, D) 
Cait 

and 

Pk) = Pkt ket) + HT (KIRK IH(K) 


Orr using the matrix inversion lemma again, 


P(k) = Plkikel) = P(kik-1)H Ck). CHCKIP(kik-1)H (k) + 
+ RCKII LHC KP (KI ko2) 
or 
P(k) = (1 = GCk)H(kK)) PC kik=1) (eyare) 
where 


By = PC ki k=l) Ho Ck) (HOCK) -PCKik=-1).H (ek) + RUK)) 


(3.19) 


From Equation (3.19) one obtains - 
Gk) = P(kikel) Hi (k) RU). 


CeCe Same) iam de) 8 ee 


Premultiplying by Sey). and using the expression found 


earlier for aaa) gives 


B{k) = Plk).{1 # Hh ROHP(Kike1)IH Ro 


[HP(kiket)H Rt + ryt = 
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1 


= P(k). CH RY + H OROHPCKik-1)H ORO. 


ee MCRD mon ae nape 
and finally, 


Bee) = PCR) .M Ch) RE) Exmauy) 


Applying this result to Equation (3.17) gives 


K(k) = KCkike-l) + GCk).Cz(k) - hx €kikel),t) 1 


(3.e1) 


If more terms were taken in the expansion of ACxCk ety ) 
than the ones considered in Equation (3.14), different esti= 
mates and covariance matrices would be obtained. Higher 


moments of the conditional densities would then be reauired. 


We are now IN Position to restart the onesteop predic= 
tion of Section III,8 and orocess a meaSurement ocurring at 


time et ° 


The initial steo in the whole process 1s at tS where it 
is assumed that an approximate GausSian distribution exists 
for the initial value of the states, x(0). Choosing the ini- 
tial estimate as the mean value of this distributions, one 


has 
x(0) = Elx(0)) 
Ele(0))] = Efx(0) = x(0)) = 0 


3V5) 





“ 


P(O) = Ele(O)e'(0)) = Var{x(0)} 


the <set of Equations) (35.7), (3.807 CSG, (3.19) and 
(3.21) is generally known as the Extended Kalmam Filter 
equations. 


D. PARAMETERS FOR THE TRACKING PROBLEM 


In order to apply the above’ relations to the passive 
tracking problem, a series of vectors and matrices first 


must be determined. 
1. System Dynamics 


Applying definition (3.6) to Equation (2.10) yields 


ie (3.22) 
i 0 Atecos BCk) = ateS(k)esin 6(k) 
0 1 Atesin b(k) At.s(k).cos 6(k) 
0 0 ! 0 
0 0 0 i 


where At = (Ct -t ) 


Breaking up Equation (2.11) into its column vectors 


gives 
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= at? nCOsS b(k) 


a (x(k)rt et) =] dat?.sin ok) 
At 
0 
~ Fat’ .8(k)esin b(k) 
g (x(k rt ty) = Sat’ .8(k).cos b(k) 
0 
at 


Definition (3.10) requires 
A (k) = 
i 


ao 83. 00 - Sat’ sin b(k) 


0 0 0 zat sCOS.b.UK) 
0 0 0 0 
0 0 0 0 


Sf 


(3.¢e5) 


(3.24) 


(S325) 





APD = CSc) 


0 0 -S At? .sin b(k) - Zat? .8(k).cos 6 (k) 
0 0 Sat °cos b(k) - dat? .s(k).sin b(k) 
0 0 0 0 
0 0 0 0 


Assuming indeoendent random forcing inputs yields 
ok) 0 
S 
Q(k) = Elwlk)wi (k)) = (3.27) 
Do 
0 (ik) 
"b 


Equations (3.8) and (3.11) require the product 


g(x(k)et rt ).Q(k) .egl (x(k) et rt ) = 
oe 6C<i Cid‘ aa ae reed" Iie 


a b d e 
b Cc 1 J 
= At. (3.28) 
d oo (k) 0 
e 0 o (kK) 


where 


a = Z At’. los 6(k) eg 2(k) + S (k) sin BK) og 2(k I 


— = At’ .sin 6(k).cos b(k) «foo Ck) - S° (kK) og 2K) 
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c= Fat “.tsin-b(k).g-Ck) ¢ s“(k).cos@b(k) .o0(k)1] 
and 


Ais = At.cos Sites cee Cie) 
S 


it) 
it 


rer a 4 2 
- 7 At.s(k).sin bCk) «0, Ck) 


me eet blk ).c-(k) 
Z Ss 


j = = At.s(k).cos b(k)ea2 Ck) 


If the expanded state variable formulation is used, 


the matrices Given below are needed 


Q(k) 1 0 
3 
| 
a°(k) = Elw?(k)wi(k)) = 1 0 (3.29) 
| 
Be ae igre 
0) 0 | of 
g+0.97 : 0 
| 
a Ge) ee pale | (3.30) 
ee —— eee oy eet 
a 2 
0 care 
1A Of 
o ran) : 
| 
| 
o°= BoB) 
yO 


3)\o) 





ee Measurements 


In the following development Equation (3.15) 1s 
aoplied to each of the measurements described in Section 


II,D to yield the appropriate linearized observation matrix. 
ae Bearing Measurement 


From Equation (2.24) one has 


(By, eek Sor Ge os, Pf © em, 0S 
7E i it af i 
aon (Cy = y. )/ Ox = x, )J 
DE i ‘i 
and then 
ee ene ete = 
Snes Ox i- DS 
=| Gothen - y,(k)1/o Gee eee 60174] 0 0 


Gs 29) 
where 
—_ 2 2 — ‘ } 2 4 ry 2 
Qa = d Ck) = (x€kikel) - x Ck) + fy(Ckikel) @- y. (k)]) 
1 
CS 5550) 
and d. (k) 1S the predicted distance of the target from buoy 
1% 
Db. Frequency Measurement 


From Equation (2.25), 
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24 = te/{Leeaas ale ort ee fcc / (c= 5. )*}. u Ss = 
Acieaet i 3 O O i: P) at 
Be ee Oe) 4 
1 Ome i Oo o& i 
where s_ is given in Equation (2.26) and 
au 
3 —_- «= =~ _o 
E84 = cos(b b+ oF s + 
+ sssin(b - b.).--b - 
alt J& 
- sesin(bd - b yo ate 
i 0& i 
Cs 549) 

and then, 

eck) = =F. (xCeik-1),t ) = 

Gat ey LS k 
= [Weert d-y 007 H(x(kik-L-x (kKJI7E By 7aA 
1 ale ae O 
(S255) 

where 


“a A 


f = f (x€kikel),-t. ) = f ECkikerl).c #7 le + 
sk i kK oO 


+ S(kik-L).cos(b(kik-1) - b Ck) 


“a 


b(k) = sretan[{yCktkeD > 97 W3DNY 


Ix(kikel) - cor] (3536) 
7 fo oS(kikel) ssin(b(kik-1) - 6 (k)) / 
1 
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( f Gk ek 
O 


A 


WD 
Tt 


§ BG, 6s ti (8558))- 


Ce Frequency Difference 


From Equation (2.26), 


mgs £.) = f/e = 2. (5. 
bo j O Ve fi 





where Ae oy 1S given in Eauation 
1 


meee) = —- (Cf. - f ) = f /e 
a O 
x(k|k-1) 


h = a..ly(kikel) = y. (k)] 
AG 1 


=|) ie) 


Measurement 


CS 4s Then, 


- a .fy(kik-1) 
J 


- f-.cos(b(kik=t) - 6 (k)) / UF (kik-1).¢) 
1 a O 


hod eC3 57) 
4 


- Yaka) 
J 


- x Ck)] 
i 


- b (k)) 
L 


h. = a elx(kikel) = x (k)) © a .lk(kikel) 

ie j j a6 
rh cos(b(kik-1) = 6 (k)) = cos(b(kik=1) 

J 
h, = s(kike1).(sin(b(kik=1) = 6 (k)) = 
1 
sin(b(ktk=1) = b (k))) 

a, = s(kikel)esin(B(kik-1) - 6 (k)) Zy 
: al 


b. (k) as in (5.56) andg as in ( 
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ra) Frequency Ratio Measurement 


From Equation (2.28), 


EN Mee ime (SS 
Bese 


Comparing with the frequency difference measure 


ment, one has 


eck) = =2 (fF 7f_) = 1/f .H_ Ck) (3.38) 
er ox ey a O di 
x(k|k-1) 
e. Time Delay Measurement 
From Equation (2.36), 
O9T 2 1/¢. Sd <-d) 
a& 13 0& i j 
where 
Bl = CoS 2G + (y - y )4 
it ie at 
and then, 
fee (k) = BCE elec= Wi/icie cen Uh h 0 0) (3.39) 
-L1y ee SIN dk 2 
x(k|k-1) 
where 
h, = (x(kikel) = x,€k)) / dik) - 
at i 


ieee ai) = Ce) oy oi) 
j j 
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Pars Vem CN) 7 ep e((y) 
i i 


Coe eon) y,(K)) / d.(k) 


Given tin Equation (3.353). 
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IV. SPECIAL CHARACTERISTICS OF KALMAN FILTERS 





Equations (3.7) and (3.8) characterize the one-step 
prediction of an Extended Kalman Filters Equations (3.18), 


(3.19) and (3.21) represent the estimation update. 


When implemented in a computing device, the timing of 
these steps 1S Shown in Fiae 5e- The actual state 18s 
represented by a broken line and the output of the filter by 
the solid linee In generals, one is dealina with vector 
processes so Figure 5 iS representative of only one com= 
ponent. Nevertheless, this. reoresentation is heloful in 


visualizing the steps involved in the estimation process. 


In a general problem the times of occurrence of meas= 
urements are not equally separated and are not Known In 
advance. At a time ty when anew set of a measurements 
becomes availabler the output of the filter 1s still the 
last estimate made close to tad unless predictions are made 


at regular time intervals independent of the spacina between 


measurements. 


A time period — 1s spent in a prediction pohase to 

p 
evaluate Equations (3.7) and (3.8), after which the best 
information about the state of the plant 1s the oredicted 


value x(keltke2). 
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If the times of occurrence of measurements were known 
in advance, the predictions could be made any time the syse 


tem was idle waiting for the new measurements and “T could 
: p 


be made zero. 


Next, the q measurements are processed using Equations 
(3.18), (5.19) and (3.21) and a time i iS spent which is 
relatively large mainly because of a matrix inversion of 
dimension (qx aq), tin the again equation, (3.19). Only after 
this time has passed is the new estimate xCce i) available. 
This estimate will be the output of the filter until] t, when 


k 


the process 1S repeated. 


Some important points in this estimation process” must 
be stressed: 

“-- the new eStimates are only available TIT + T 
seconds after the arrival of new measurements. 

“-- a major part of the time period of duration i 1S 
spent in the gain equation due to the required matrix inver= 
SION. 

“--- from the end of the computation of an estimate 
until the arrival of a new set of measurements the computina 
equipment 1s idle and thus free to execute other chores. 

“-- the one-step prediction is based on a linearization 
of the plant dynamics over an extended period, for example, 
from t to t_. 

k-1 k 
“-- the measurement equations are based on a Jineariza= 


tion about the predicted values. 
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A. PARTITIONING OF MEASUREMENTS 


When some of the meaSurements come from statistically 
independent sources, steps can be taken to alleviate orob- 
lems which arise. Considering the special case when all the 
measurement noise components are independents, the covariance 
matrix Rk) is diagonal and the following develonment 


applies. 


1. The Linear Case 


Assume a dynamic system with linear observations of 


the form 
Zak ) = Hk) .x(k) + vk) (4.1) 
where H(k) 15 not a function of the states. 


The estimation eauations of Chacter III become the 


Standard Kalman Filter equations 


mae PCki k~1)H Ck) (HOEK) PCK 'k-1)H Ck) + Rk) IL (ie) 
Beers = x(kik~1) + GOk)iz(k) = H(k).x(kik-1)} (4.3) 
P(k) = (1 = GCK)HCKIIPCKIk=1) (4.4) 


The measurements are assumed to occur simultaneously 
and, as shown for the first time in (9], for this linear 


case, the results obtained by processing all q measurements 
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at once using these equations are the same as the final 


results obtained by processing one measurement at a time. 


In the latter approach, the result of processing one 
measurement component is used in the following computation 


to process the next measurement component. 


A general way to show that the two approaches) yteld 


the same estimate after al] measurement components have been 


processed is given in the following paragraohs. 


In shomt it is to be proved that the estamates qaaG®) 
and P(k) obtained by the use of Equations (4.2), (4.5) and 
(4.4), when all measurements are grouped in a vector Zz, are 
the same as those obtained by the use of the following 


iterative equattons aq times, once for each components: 


Eee= Pe GH 7 HP OA? eR (4.5) 
i i-l aft iil i i 
= $ = = 
mee 5) Og 2g 
= (I = G.H.Jx. _+ G2. (4.6) 
ii -i-l ii 
PB = {I -GHI.P ee) 
i ii -l 


1 = 1, cr eee ff Q 


where 


P = P(kik=l), ee eee) 
fe) —o = 
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—~q 

2 H V 
il 1 1 
z. Ho V9 

20k) =}. = ; x(k) + . 

NN. 
2 H v 
q q q 


If Equation (4.6) is applied a times,s the following 
relationship between x(k) (or sf ) and x(kikel) (or x ) 
— = = 


results 


Sawer = Caroli toeeG Wy dimen (he=) GH Ix ch) esioue 
a q q Ga Ea g— 1) aR 


+ G z2 + (f = GHIG 2 + wee + 
qq @q g-l-q-i 


eo Ghia (l o Cudted cool = Eh A, 


(4.8) 


and recursive use of Equation (4.7) gives 


et{k) = Hi = G HIE = Ga tgJeee Ul = GH IP Ckiked) 


(4.9) 
Comparing (4.8) to (4.3) and (4.9) to (4.4) one sees 


that the assertion can he proved by showina the validity of 


the Fevat ions 
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[Ll = GC(kJHCkK)) = bf =6. JA) = Giese 
qq iS 


q-l ql 
(4.10) 
| 
(n x 1) eCopt) (nr x 1) 
ae 3 
mmee=i(] - GH J... fl = G.H.JG. !...! {1 = GH JG ! G 
) Sf 2201; | qq ql q 
| 
ce ee 
Or, in & compressed way, 
(n x Caqe1)) eC racan la) 
| 
Gtk) = |({l-GH1cGe | G (4,12) 
q q | q 
and 
Be CK HICK) IBS (i = Gan} (ie 6.47] (4.13) 


* 
where G is the gain obtained from Equation (3.19) if ane 
has a q-1 dimensional] measurement vectors and H* and R*are 


equivalent matricesr as shown below 


aie 
H(k) = |-- - (4,14) 
lH 
q 
R™ 1 9 
| 
R(k) = Je ee bh - (4.15) 
| 
‘ae ie 
1 q 
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G* = PCkik-1)H* CH*PCKik<t)H* + R* A (4.16) 


Let us expand Equation (4.2) and try to manipulate 


it into the form of Eauations (4.10) = (4.15). We have 


I it -. 
i i 
T | At 
r | 
ee =e = | = = 
; 
eee i eee { eee j eee 
} ae 
: on 
Hapoee b eee eee kh PEE a 
| ne fae de a 
{ | t 
or 
i 
a ! 25 
[HPH' +R} = 235 ale sos (4.17) 
T 
23 | “4 
where 


A, = te trey ees 
= wuxtpyl 
oe at 
er go. fT 
Se se = a, 
a, =H PH +R 
4 qq q 


Next define the scalar ¢c as 


c= a,- afta 
Seer 2 


ve 





Thus, 


ec =H PH!’ + R «~ H PH! cCH* PH ls R*p4y* PHI 
qq q q q 


Using the definition of G in Equation (4.16) gives 


e 


ec =H PH’ +R «~ H G*HtPeH! = 
dua q q q 


= H [I - G*H*)] PH! + R (4,18) 
q q q 


It is shown in [10) that 


I i | 
A B 
| = ieee ee = 
I = I 
| 
T | ay I 
53 3} 
= =i 
where 
—~ awl -l T -~L_-l 
BF A ft A, ana,A cs (4.19a) 
_.,l. -l 
b, = A, a¢ (4.19b) 
T T.-l -l 
= - 4. 
b, are (4.19¢) 
~ Ll 
b = ¢ C4719 ae 


Using these relationships, 


G = PH (HPH + RP? = 


i 
= eT, = 

ES BM at 
-= LJ ie > 

| 

| 

r 
3 1 8, 


(ie 





*T 


crPH te + pH'ph! +: pH*ly + PHIb? = 
q-3 oe q 4 


1 


(C ; 0) (4.20) 


Applying Equations (4.19a) = (4.19d) to (4.20), 


* = = 
(H™PH SS + RY tH eHTSly pH l= 
q q 


GT H*PH PH ele (1 = G*H*) PHT cls 
q 


q 
Ak Jb ee, 
= TMemcoH | Pen? = Gon )eHi ee Re ee 
q q q q 
Let P* = [I = G* H* JP be the estimation error 


covariance matrix after the processing of a=! measurement 


components, as seen from Equation (3.18). Then, 


Oe P Hope p“H! + e Fle c (4.21) 
7 4 qed q q 


1S the gain when processing the qth measurement. 
Ne also have 


* = 
C= PH CH“PH? + Rey + PH*T (WH *PH*T + R*PlH* PHTH PH*T 
qq 


x #e = = 
(Hic Eine jen oor Comp ler pa lege etre lee 
qq 
* = = 
Sco GG Heed HoG co. = ereame- —! 
qq qq 


Since c is a scalars one can write 


E = = nll = Gabel eHuc nce = Cem- Drie” 
qq q 


and, therefore, 
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@ = tS carene (4.22) 
: q q 


and from Clara 
* { 
G=j)(1 -GHIG {| G (4.2e5) 
% qq ! q 


It is also the case that 


x. x 
GH = Ec - GH JG | H 
qq qd 


H 
To 


=e, [ee cmc H+ cup 
qq qq 
SO 


m= GHleet - GH - (m= CM GG -H* 
qq qq 


and, thus 


[-GH = (1 - GHC - a (4.24) 


Equations (4.23) and (4.24) are the same as Equa- 


tions (4.12) and (4.13), proving the initial proposition. 


Computationally,s, one can replace a {axq) matrix 


\ 
inversion anda (nxq) x (qxq) matrix multiplication by only 


q scalar divisions. For an indication of how much computing 


time can be saved, consider that a matrix multiplication 


- scalar multtolications and addi- 


3 


requires approximately na 
tions, and that a matrix tinversion requires at least q~ mul@ 


tipnlications and additions, not counting logic and indexing 





time nor considering storage requirements. 


The partitioning of the measurements then 
corresponds to trading at least q*(n + q) multiplications 
and additions for q divisions. For a machine with average 
multiplication time of 10 microseconds, addition time of e 
microseconds and division time of 12 microseconds, if n= and 
q are equal to 4 a saving of about 1.5 ms is obtaineds if n 
and q are of the order of 20, this saving could be of the 


order of 200 ms at each measurement time. 


Besides this comouting time savinge, the processing 
of one measurement at a time provides partial estimates that 
can be used as improvements over the predicted values, as 
Shown in Fiqure 6. Also, the final estimate is obtained 
sooner, although with the same accuracy as when the measure= 


ments are all processed simultaneously. 


ce The Nonlinear Case 


We assume the same system as before but with non- 


linear observations of the form 


X 


z(k) = (xCk)rt toa CK’) 
By linearizing this equation as shown in Chapter 
Ill, using the definition 


A) “a 
H(k) => — cee lel eg ’ 
GE) See Ina kei oad 


\>< 


the proof of the last section 3s still valid for the 
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a. —-| X (k) 





ate seit +T +T time 
k ate Pp me p e 
q predicted value end of processing 
measurements available; start all q measurements; 
available. of processing all estimated value 
q measurements. available. 





SUC TC= SERIE SEE 
A - time 
start of end of processing all q 
predicted Walu 


roce 
processing measurements; final estimated 
second meas. ; 
available; ee. rE value available. 
of processing 
first measurement. 
Figure 6 : Partitioning of measurements = 


Linear case. 


a 





linearized equations. 


If, however, the partial estimates obtained by pro- 
cessing one measurement Chopefully a better approximation to 
x(k) than CK ety) results) are used to calculate the A 
vector for the processing of the next measurement, better 
accuracy should be expected than when al/| it vectors DIG 
obtained using the relatively crude predicted value. This 


effect should be particularly significant when the time 


between each group of measurements is relatively large. 


ThusSe partitioning the measurements in an Extended 
Kalman Filter should provide improvement in estimation accu- 
racy in addition to the advantages pointed out in the previ-e- 
ous sectione This imorovement can be visualized as shown in 


Figure 7. 


The expected improvement in the estimation process, 
brought by partitionina the measurements in Extended Kalman 


Filters, can be explained in terms of the conditional densi- 


ties. 
Consider two independent simultaneous measurements, 
h = 
Zz, Ck) 1 620K) vi Ck) 
z(k) = = 
Zz. («k) h (x(k)) + vi Ck) 
2 We 2 _/ 


The conditional aensity after processing the two 


measurements iss from Equation (3.1), 
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T mele" we ee we 
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processing the 
q measurements 
together. 


=" 
(k) | 






vec ie 6 6 saat s 


_ ee @@8 7°33 2S & ee S| SS = a ed 


-T 
Se 
end fof end of processing 


PEOCessinge sc q measurements 
measurements alltogether. 
one at a time. 
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Figure 7 : Partitioning of measurements - 
Nonlinear case. 


is 





pcieuyzZ*) = pix(k) 12, (kez, deo = 


plz, Ck) e25(KIIX(KI) POCKIZ DY /plz, (kez (KID = 


t ry -l 
(p29 CkI EZ (KI ex(KI) / p(z, Ck) iz, (k) rz" re 


fle, OD HCO) CLA We CED Ae 


But, from (3.1) 
p(x(k)iz,(kIeZ) = plzz (kIEC II CK(KIIZ YH 7 
pz, (ki Ze) 
also, 
P(zo(k)iz, CkIexCkI) = Oz, Ckdix€k)) 


because of the independence assumption. 


And so, 
o(x(kIIZE) = tol2 (kIix(k)) 7 lz (kIIz (kIVZ J, 


p(x(k) $2 (ke) (4.25) 


The terms on the right hand side of (4.25) are: 


p(x(kitz, (k),Z4) “-- has all the information about the 
States from processing up to the first measurement component 


zi Ck). 


Seek) px CK) mea S Obtained  Grrectityemtrom the 


Z 
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measurement equation and doesn't depend on the first mease- 


urement. 


aba Usb trs Oso --- is obtained from the measurement 
equation but using all the results given by the processing 


of the first measurement component. 


This last term indicates that information is lost 
if, ina Extended Kalman Filter, the linearization for pro- 
cessing a measurement iS based on the predicted value 
Igexweao Of on the value estimated after processing the 


oreceding measurement component. 


Se The Tracking Problem 


In the tracking problem here addressed, the measure= 
ments naturally occur at different times. Bearing indica= 
tions are basically continuous or, if preprocessed, succes= 
sive measurements are available within short intervals of 
iigess Frequency indications must be obtained from digital 
processing of records of considerable time length. Time 
delay indications are also obtained from diaqital processing 


but with longer time intervals. 


In earlier work, [1], [2] and (3), these measure= 
ments are forced to occur simultaneously and are processed 


together. 


In the previous sections of this chapter it was 


shown that, if the measurements are all available at the 
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Same time, they should be processed separately for better 
computing efficiency ands with nonlinear measurements, for 


possibly better estimation accuracy. 


Since one 1S now free from having to process” the 
measurements together, and encouraged to do Sor, why not use 
the idle time of the computing equipment and process” the 


measurements at different times? 


If one does this the decrease in computation time 
obtained previously will be diminished because of the extra 
computation time required to do many predictions, one for 
each time a measurement 1S to be processed. Nevertheless, 


the advantages are potentially of great value, that is, 


-=-- a measurement can be processed as soon as it 1s 


availabler, with no waiting for other measurements, 


\“--- measurements that occur more often can de pro- 
cessed more times than measurements ocurring less free 


quently. 


“-- the output of the filter is updated frequently. 


m-=- the one=step prediction 1s based on linearizae 
tions through much smaller time=steps than before, hopefully 


with improved accuracy and fewer divergence problems. 


Figure 8 gives an idea of the Improvement one 
expects to obtain by (1) processing simultaneous measure= 


ments separately and (2) processing the measurements as they 
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naturally occur, when a Extended Kalman Filter is being 


used. 


B. GRAPHICAL INTERPRETATION 


Benet ail understanding of the actions taken by a Kalman 
Filter during the erocessing of a meaSurement can  0be 
obtained through the graphical interpretation and visualiza= 


tion of these actions in a twor=dimensional problem. 


e 


From the conclusions of Section I[V,A, it is advantageous 
to procesS any group of measurements with statistically 
independent noise components Separately. In this section, 
thus, the processing of only one meaSurement at a time is 


considered. 


ie, The Linear Measurement Case 


We will consiaer a dynamic system with two state 


variablesrs anda linear measurement of the form 
z2(k) = H(k) x(k) + vek) = 
= hy Ck) ex, Ck) + h, Ck ox, Ck) + v(k) 


with v(k) @ zero=mean,s, discrete white Gaussian noise process 


with Elv-(k)) = r2(k). 


At time th a prediction 1S made from previous esti-= 
mates and knowledge of the olant and system noise Charac= 


teristics. Equations (4.3) and (4.4) are used for a linear 


&3 











measurements together. 


eee 
-——_—_—_— 
| ae a 
= = o_o 
| ’ _—- J ae ae 
k) ae ' | 
x(k) ee processing simultaneous ft eas 
ea / measurements separately. 
A —-—- J 
X (k) SL aa ee 
: 
, st processing simultaneous 
| 
| 


ee 
re Stil 






processing heasurements 
as they naturally occpr. 


“> | 
t] 


Figure 8 : Comparision of processing policies. 
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system, (3.7) and (3.8) for a general nonlinear system. 


The Kalman equations (4.5), (4.6) and (4.7) are now 
applied to correct the predicted value according to the new 


measurement. 


Let's represent the Quantities and relations 

involved as shown in Figures 9 and 10% where ei and e area 
: | Z 

pair of orthonormal basis vectors and a general vector is 


represented by 


The following conventions and simplifications apply 


rs 
x, Ok) x* 
xx 5 x(k) = = 
x(k) x * 
2 2 | 
x 
1. 
x = x(k) = 
2 
a 
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actuals state 


predictad state 





Figure 9 : Graphical interpretation - Measuremeat 


lines and gradient. 
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Figure 10 : Graphical interpretation - Error ellipse. 
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1 
g = G(k) = = 
g, Ck) 9. 
ee Peel / 2rd 21/2 
fhp = Cho hy’*= Cho + Woh 
mee” Ze) 2 
fof = (9, + 95} 


= Sh 
1 arctan Ch, i? 


= / 
B arctan tg, 9,) 


Let's define the function 
@ (x) = Aix, + hy Xx, (4.26) 


and call “measurement line” the locus of points where $ Is 


constant. Some of these lines are shown in Fig 9, 


Using this definition one sees that the previously 


defined h 1s the gradient of 9 with respect to xs or oe , 


which in this linear measurement case 1s not a function of 


the states. 


The prediction error covariance matrix,y 


Pik kel jes Pl = 
o! 
a2 Lp 
1S represented by an error elliose whose direction (6) #4=x®and 


axis dimensions (a = major axis 6 = minor) are given by, 
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(2) 


Ms = 2 § / ¢ = ’ e.6 oe. = < Pa < G.ef7 

tan 28 Po (5 256 T 8 T ( ) 
= ' 5 2 : + es fe ine © A226 

= (Py Doe aD og 
my ¥ + : / ~ : in 28 4.29 

b (Po, P45? e 2a $1 ( ) 


ae the State Estimate 


Equation (4.6) can be represented in the form 


| <> 


= x' + g.fresidual] (4,30) 


where 


(residual] = z(k) = H(k).xCkikel) = 
=o (x*) =o Cx") + vk) CRESS 


and (4.50) means that a vector correction of magnitude 
l\g{.{residual} is made to the predicted state in the direc- 
tion of the vector a which has an angle g as shown in Fig-= 


ure ll. 


Equation (4.5) gives, as shown in Appendix A, 


on ee De at 
2 De i) a 

2 2 2 

p' ho + 2@pt hh + po' fh +r 

Sy an Hl at bay fe: 22 2 
q = = (4.52) 

Ye ho) sane 
95 aie LD 

ony - he? pee 


t 
1% * Sy ALR + Poh * SF 
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o(x) = > (x*) 


d(x) =\9$(x') 
XF¥e—------- —Ke— 
| 
3} 772 EX 
n ys" a B 
ve 
) — 
pesal 
| | 
| | 
| | 
e, ny 
| | 
! 
2 >. se 
a il il 


Figure ll : Graphical interpretation - Gain vector. 
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(1) The Direction of Correction 


The direction of g is seen to be tndepen- 
dent of the measurement noises orr given h and P(kik=-1) the 


direction of correction is already determined by 


= q ¢ q q 
B arctan UCp', hy + P55h,) / (pi,4, + 5 hI] 


(4.33) 


From Equations (4.27) = (4.29) one obtains 


t bend a = e : (4.34 
Pio : (a b).sin 28 ) 
o' 2st ((a +b) + Ca - b).cos 28] (4.35) 

lal 2 

so si - ~ bo). 9 Ae 
P59 5 [Ca + b) (a b).cos 28] ( ) 


Apolying these equations and the definition 
of a to Se oe) gives, 
a7 = 1).sineo + (Ca/b + 1) = (a/b = 1).cose6).tana 


Se ie 
Ca/b + 1) + €Ca/ob =- 1).coseds + (a/b e-1).sinedd.tana 


(4.57) 
which can be reduced to 
(a/b).tango.e([tangetang + 1) + [tana = tané] 
a 8 te 
(a/b). (taneetang + 1] = tang. (tana = tanég] 
(4.38) 
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This equation shows that the direction of 
the correction generated by the Kalman Filters tn this two- 
state problem, is a function only of the direction (aq) of 
the gradient (h) and of the alignment (6) and the ratio 


(a/b) of the prediction error ellipse. 


Three special cases of interest occur when 


(i) a/b = 1, (ii) b = 0, and (iit) G6=a . 


Case (i) = If a/b = 1 the error "ellipse" 
is actually a circle meaning that the prediction has the 
same uncertainty In any direction. No preferable direction 
pre-exists and the correction 1s in the direction of the 


gradient of the measurement functions as can be seen from 


(4.57) 

tan 8 = (0 + e.tana) / (2 + 0) = tana 
or 

B=a+t+nwmT 

Case (ii) - If b = 0 the error ellipse 1s a 

line meaning that the prediction 1s exact along the minor 
axis and that only corrections along the mafor axis. are 
necessary. As can be anticipated, the direction of correc 


tion will be the direction of the major axisSs as can be seen 


by taking the limit as b>-0 of Eq. (4.30), which yields 


tan 8 = tané 
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or 
8B = O6+¢ Ny 


For any other intermediate value of a/b between 1 and 00 ’ 


the value of 8 will be between qa and 6 (+ no. 


Case (i111) = When 6 = q the gradient is 
colinear with the major axis of the prediction error ellipse 
andre from (4.38), 


tang = [(a/b) .tan6.(tano+ 1) + 0) J/ [(a/b). (tare + 1) = Q) 


or 


tang = tanée and 8 => @+ An = atonan 


Thus the direction of the correction is 
along the major axis of the ellipse which 18 colinear with 


the gradient of the measurement function. 


(2) The Amount of Correction 


The magnitude rele the correction 
({residual|{.ig{) can be better seen if one decomposes a into 
components m and ne along the gradient and the tangent to 


the measurement functions as shown in Fig. 1l. 
ms q eCOS a + g SING (4.39) 


1 


and from the definition of az 
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Go's’ o Seah. 7 {hile Sima = haji (4.40) 


Usine Equations (4.32), (4.39) and (4.40), 


5! h2 + p’ hho + pop’ hh + O° h 2 i 

mee, Ll 1 2 py ae 
La aaa h2 a h 
eat Ne Soe Ee ou 


= ay oe (4.41) 


where 


yee s/s (pt he + apt hh + pt he) (4.42) 


i oma 12 thee 22-2 

Consider some special casese How would the 
new estimate be computed if a noise-free measurement were 
obtained and the Kalman filter were aware of this fact, 


1eCes ro 0? 


Note from Equations (4.41) and (4.42) that 


when the measurement is noise-free, r = 0, Y = 0 and 
msi / {tht , @ constant. 


No matter what P(kik=-1) is, the projection 
of the correction into the gradient is a constant if r = 0, 
or, the tip of the correction vector follows a line normal 


to the gradient when P(kik=-1l) ts varied. 


Figure 12 shows the correction made to the 


predicted state when b = 0. In Appendix C it is shown that 
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for the present case of a jinear measurement function the 
vector OA, in the direction of h and length [residual] /jh¢ 


has its end point A on the line 


¢(x) = o(x') + residual = g(x") + go(x*) + v = > (x') 
or 

$(x) = o(x*) + v = z 
which, for r = 0 (i.e. v = 0) passes through the point x = 
x*%, 


Since the correction has the direction of 
the prediction error ellipse (line) and has OA as its pro- 
jection into the vector hs, it can be seen from Figure tie 
that in this case of b = 0 and r = 0 the estimate will coin- 


cide with the true state. 


Fiaqure 13 shows the correction for a gene 
eral error ellipse but still with r = 0. It 18 seen that the 
projection into h 1s the same as before. The direction of 


correction 8 is Somewhere between a and6. 


A general correction for a noisy measure=> 
ment is shown in Figure 14, with different scaling than the 
previous figures. The direction of correction g is first 
determined and suppose it is as shown. The projection of the 
correction into h would be OA if r were zeros for r £0 the 
reduction factor is applied and suppose the projection is 


vector OB as shown in the figure. Taking a normal to h- from 
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(x) = z 


error Je iin sew leiency 
(x) =\$(x') 
perfect estimation 
Zo 
- J | 
2 —. g-.-residual 
7 


\ 


length = residual/|{h|{ 


- 





Figure 12 : Linear measurement - Perfect estimation 
EOnso—Umand: <—O. 
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o(x) =z 
¢(x) =\¢(x') 





Peres tduae 


x 
ya residual / |h| 
A 
x, "| - B 
x! 


Figure 13 Linear measurement - r=0. 
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point Br, the point C where it cuts the direction of correc 
tion determines the final correction OC and the new esti- 


mate. 


b. The Estimation Error Covariance Matrix 


The effects on the estimation error covariance 
matrix can be seen in the following waye From Equation 


(G.7), 


P(k) = (1 = GCK)HCK)IPCKik=1) 


or 
i = GWA SY es 
al Pio ar 1 PT 2 an Pi2 
9 ne eee eo p' p' 
12 22 21 22 12 22 


Using the values of or and a5 from (4.32) one 


finds that 


= Z g § 2 t Z 
P4117 Ch, (p45 P45 Py'5 ) + oe r ) /e (4.43a) 
a 2 2 
= (hoh Co! ~ : : + : / (4.43b) 
P12 iO opmiigs lain, ipey | Me es 
a i Q ‘ 12 ’ 2 : 
pop = (h (p' p - 9" ) + 9° r4) Se (4.43c) 
2. i vit 22 i? ZZ 
where 
= Z. Z 2 
oman © oan do ee '" Hho + } (4.43qa) 
Sean [> 1 2° “oo 5 
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d(x) = o(x*) 


> (x) O(x*) +v = z 


(x) = O(x') 






aoe , _-residual /|h| 
g-residua 


vA 
vA 
A Bor: a 
x, A 
m.resNdual 

t 
=? 
=9 


Figure 14 : Linear measurement - General correction. 


Ly, 





And the estimation error ellipse has the orien 


tation 


tan 26 = CP 1 o/ (Pi, - Poo) = 





2 2 
e = ¢ 9 + ! 
eal ly 7 inp, Ya anal uae 
Ro ap Ci a es ee 
(hy 2) 8Pyy = Py Poo? Aa Jog Ul 


The limiting situations occur in the special 


cases when (1) r +o and (ii) r = 0. 


Case (1) = When r is relatively larger little 
information is provided to the filterr the estimate is 
essentially the same as the predictions and the estimation 
error ellipse is essentially the same as the prediction 
error ellipse. This can be seen from Equations (4.43a) ° 


(4.430), as Pr ->} © 


Case (ii) = When r = 0, all errors in the pred- 
iction, in the direction of the qradient of the measurement 
function, are corrected and the estimation error ellipse 
becomes a line colinear with the measurement line. As shown 
below, this result is obtained from Equations (4.43) and 


(4.44) when r is made equal zero. 


From (4.43a) - (4.43d), p .p is easily shown 
ale i 22 


2 
12 


singular. From (4.44), 


to be equal to p when r = 0, indicating that Pk) hecomes 


(ho - nh“) = 


tan es = chi h, / 1 9 
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2th s/h.) / C1 = Ch sh 24) & 
Fe hig il 


Dene? (Ole Gago) 2 bap ae 


From analysis of the individual] terms of (4.43a) 


= (4.43d), this result 1s to be interpreted as 


20 = 2a + 7 


5 = a ¢ m/e 


When the minor axis of the prediction ellipse is 
12 ae q q 
AID 2 Va 


measurement 1s noiseless, then P(k) = 0 and comolete cere 


already zero, as in Figure 1iecyr i.@e, , and the 


tainty about the system state 1s obtained. 


These two last oaraaraphs Show a very interest- 
ing situation: if two oerfect measurements are made on a 
twoestate plant and the filter is aware of this (j.e. r = 0 
is used in the Kalman equations), then an exact estimate is 
Beeoned ~-~ thie estimation error covariance matrix col- 


lapses into the zero matrix. 


For a general situation with normal measurement 
noise values, the estimation error ellipse is rotated toward 
the measurement line, away from the gradients with 6 occue 


PYING a position between 6 and qo te TI. 
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The general expressions for the lengths of the 
major and the minor axis of the estimation error ellipse can 
be obtained by the application of Equations (4.28) and 


(4.29) to (4.43a) = (4.43c). 


ce The Important Points 
Summarizing the principal points seen up to now, 


“--= the direction of correction in the estimate 1s 
defined at the end of the prediction phase and doesn't 


depend on the measurements. 


--- if the measurement is noiseless and r= 0 is used 
in. the Kalman filter equations, the projection of the 
correction on the gradient of the measurement function is 
(residual / |{h]),s independent of P(kik-1). The estimation 
error covariance matrix becomes singular indicating, in a 
two-dimensional problem, that the estimation error ellipse 
becomes a line **= a point when the prediction error ellipse 


is already a line. 


“-- if the measurement is noiSys the projection of the 
correetion on the gradient is multiplied by the reducing 
meactor (1 / (i + Toes The error ellipse is rotated toward 


the measurement lines, away from the gradient of the measuree 


ment function. 


e-- if the measurement 1S too noiSys the reducina face 


tor tends to zero and the estimation is the same as the 
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prediction. 


ae The Nonlinear Measurement Case 


Assume the same twordimensional system as before but 


with nonlinear observations of the form 


z(k) = h(x(k)et 2 + v(k) 


ae The Extended Kalman Filter Approach 


The prediction phase is the same as before but 
the estimation equations are now (3.18), (3.19) and (3.21) 


for an Extended Kalman Filter. 
Let's. again define a measurement function by 
g(x) = h(x) 


and redefine h as the gradient of 6 with respect to x. 


Now the components of h are no longer constants 


but depend on xe At the predicted point these components 


will be 
h = h(x(kik-l1))  , h =% A(x k}kel)) 
oe a dX er 
Equations (4.27) through (4.44) are all still 
valid. 


Let's consider two tyoical measurement functions 


and observe granohically how the new estimate is obtained. 
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Case (1) = assume the measurement function is of 


the form 


o(x) = arctan (Cx, = ai? / Cx, = a, )) (4.45) 


The measurement lines will be as shown in Figure 15. 


Suppose now a noiseless measurement is obtained 
and is to be wrocessed. Fig. 16 shows the situation when the 
prediction error ellipse is a line (ob = 0). In Appendix C it 
is Shown that for this special measurement function the vec 
tor OA of length Cresidual] / fh} will end somewhere before a 
point 8B on the curve g(x) = 2. Supposing its length is as 
shownr, 1t can be seen that an imperfect estimate is 
obtained. Nevertheless the estimation error covariance 
matrixs given also by Equations (4.43), will collapse into 
the zero matrix indicating, falsely,r that a perfect estimate 


was obtained. 


Figure !7 shows a typical situation for a geno 
eral prediction error and a noise-free measurement. The tio 
of the correction vector can be seen to follow the line A = 
A' parallel to the line (x) = g(x") instead of being along 
the line (x) = go €x*). AS seen in the previous section, 
the estimation error ellipse will erroneously be determined 
as a line also along A = A', indicating falsely the direce 


tion of the true state. 


Case (ii) = assume the measurement function 1s 


of the form 
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Figure 15 
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Nonlinear measurement lines I. 
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Figure 16 : Nonlinear measurement I - Imperfect 
correction for b=0 and r=0. 
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Figure 17 
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ye? 

(x) =| ¢ Yt V7 (4.46) 

pr x? - xa oT 7) a) e 
The measurement lines are as Shown in Figure 18. This’ same 
figure shows the new estimate obtained when the measurement 
is noiseless (r = 0) and the prediction error ellipse is a 
line (b = 0). The same observations of Case (i) are valid. 
As shown in Appendix C the point A 1s on the curve g(x) = 2z 


for this special measurement function. 


These figures have their values mostly because 
they shows, in a simple way, the large errors that can be 
made in the processing of nonlinear measurements by an 
Extended Kalman Filter. Two types of errors are generated: 
one in the determination of the new state estimates the 
other is seen in the fact that the estimation error covari= 
ance matrix P becomes a poor approximation to the true error 
covariance, making the filter non-optimal. These difficul= 
ties are a function of the nonlinearities involved and, 
apparently, are more sianificant when the predicted values 
are poor and when the measurement noise is assumed to be 


small when evaluating the extended Kalman filter equations. 


b. The Iterative Approach 


A way of correcting, in part, these deficiencies 
is the adoption of iterative procedures. The first improve- 
ments can be obtained ifs after the Extended Kalman Filter 
has been applied to correct the prediction, the filtering 


process 1S repeated but with the linearizations made about 
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o (x) = o (x*) = "7 


¢(x) = $(x') My , g. residual 


ec dale s|iay| 





Figure 18 : Nonlinear measurement II - Imperfect 
correction for b=0 and r=0. 
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the just obtained estimate, hopefully a better approximation 


to the true state than the crude prediction. 


Analytically the process consists of replacing 
the estimation equations (3.18), (3.19) and (3.21) by the 


iterative equations 


x = kCkikel) + G Cz€k) = hx Ckikel))) (4.47) 
i+] uu =o 
G, = PCkIk-1)H, (HPC kik-1)H + RY (4.48) 
a aL a 
A. — h(x. ) ge ; = 0, 1, Ce eee # f 
- 2 9x fl ll 
where 
x = x(kik-1)  , | x = x(k) 
=) — —f —_— 
P(k) = {I - GH JPCkike1) (4.49) 


For our twosdimensional problem Equation (4.47) 


can be expressed as 


A 


x, = x' +g .residual] (4.50) 
liye dt 


where 


residual = z= h(€x') = (xr) = ox") + vk) fo ey), 


Adding and subtracting (x ) to the residual, 


L 





or 


¢ = & = 4 ze ’ 
eat 9, lz o (x50) 9, - lolx) ox) (4 5¢c) 
These equations can be applied until there is no 
significant difference between consecutive estimates, and 2 


or 5 iterations are normally enough. 


The effect of these equations can be easily 
visualized for the simple case of b= 0 and r = Q. The 
sequence of Figures 19 = 21 show this effect for one of the 
special measurement functions already studied. In Figure 19 
the first estimate is obtained as done previously in Figure 
16. The vector QA, which is smaller than OB as shown in 
Appendix C, represents the projection of the correction = on 
the gradient hoe A normal to sat from point A cuts the pred= 


iction error line (b = 0) at the first estimate or 


In Figure 20 the terms of Equation (4.5e) are 


shown for 1 = 1. =j%PTwo correction terms are applied to the 


prediction x', both in the direction of 9, but with dif- 


ferent residuals. The direction of g_ is again along the 


tes 
prediction error line (b = 0). 


The gradient of the measurement function at the 


first estimate is h The projection of the first correc 


lee 
tion (g,fz = d(x 1) into h, is smaller than CE, say CD. 


Taking the normal at D one gets the first correction term 
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Figure 19 : Iterative process - First estimate. 
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CH. 


The projection of the second correction term, 

on f é(x') - V1 JAD AU taliete: me is smaller than CG, say CF. 
~- aes — 

The normal from point F determines the second correction 


term (IC). 


Adding the three vectors x', CH and IC one _ has 
the second estimate, shown in Figure els and the next itera= 


tion can be processed. 


If in Equation (4.47) one modifies the residual 


by 


mau) = h(x(kik-1))) = 20k)Be (h(x + 
1 


+ H .(x(kikel) = x.) + 2.) = 
L a a 


= 2(k) = h(X%.) © Ho fk(kikel) = &] 
a is 1 aa aoe lk 


- 
nape 5 


then Equation (4.47) becomes 


ps2 


= X(kikel) + Gitz(k) = he 


Fes 20 oR USGheih SBM 


(4.53) 
This variation in the residual is very small to 
affect considerably the estimates but in (11) it 1s Shown 
that its use produces a better agreement between the calcue 
lated estimation error covariance matrix and the true error 


covariance, when processing a noisy measurement. 
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Figure 20 : Iterative process - Correction terms. 
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5 Summary 


The direct application of these graphical interpre- 
tations to a general problem is impractical. Nevertheless 
the basic points observed here can be generalized for Kalman 


Fylterss: 


—_ 


“--- the “direction” of correction in the estimation 
is defined by the predicted quantities and doesn't depend on 


the measurements. 


“-- the amount of PORT ection in an estimations along 
the direction defineds 1S the result of a weighting process 
between the uncertainty in the measurement and the confi- 
dence in the predicted values. This correction is a function 
of the factor (1 / (1 Sa) where, in Qqenerals for a n- 


dimensional system 


n n 
we tk ete (4.51) 


i-l j=l ij i j 

“-- the errors in the estimate from the processing 

oO f a nonlinear measurement can be larae ands most important, 
the estimation error covariance matrix may become a- very 
poor approximation of the true error covariance, forcing the 
filter to tnaccurately weight future measurements. This 
problem depends on the nonlinearities involved and seems to 
be worse when the predicted values are poor and/or the meas- 
urements are considered by the filter to have relatively low 


NO1S@. 





“-- the errors indicated above may be reduced with 
the use of the iterative equations (4.48), (4.49) and (4.47) 
or (4.55), instead of the standard approach given by the 
Equations (3.18),(03.19) and (3.21), at the expense of 


greater computing effort. 





Ve. COMPUTER SIMULATION 


A. GENERAL CONSIDERATIONS 


A computer simulation was implemented in a Unix Operat- 
ing System running on a POP 11/50 with the FP 118 floating 
point processor, located at room 506 of Spanagel Hall, NPS. 
The whole program is interactive and is expected to be 
self-explanatory to a user. Its block diagrams and coding 


are presented in Appendix E. 


The araphical outputs were generated on a Tektronix 
4GO14=-} terminal. The presentation of the results of the 
various simulations is simolified by the adoption of the 


following conventions. 


The metric system was used in all calculations and out 
puts. The axes are dimensioned in kilometers. Heading and 
bearing angles, expressed in radians, are measured counter- 


clockwise from the positive X-axis. 


Figure 22 shows a typical plot. Points of the true track 
are represented by the sysmbol "+" and interconnected by a 
solid line. Buoys are represented by a distinctive boxed 
symbol. The portion of the true track that is in the range 


of a buoy 18 bounded, if necessary, by dash=dotted lines. 
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The signals received by the sonobuoys are assumed to be 
sent to a ship or aircraft where they ar® processed by spe- 
cialized equipment to generate certain measurements. The 
major characteristics of these measurements are printed on 


top of each figure and obey a code whose aeneral form is 
L Ni Ne NS § 


where 
Ll =- a capital letter indicating a type of measurement. 
B Stands for bearings FF for frequency? D for frequency 


difference; R for frequency ratios and T for time delay. 


Ni = an optional number indicatingr when necessary, the 
buoy that generates the measurement. If this number has two 
digits it identifies two buoys, thus, T le indicates a time 


delay measurement between buoys 1 and e. 


Ne - anumber indicating how many seconds) after the 
start of the problem this measurement will be first 


obtained. It is also optional. 


N3 = a number indicating the period between consecutive 


measurements. 


S =* a number indicating the standard deviation of the 
simulated zero-mean Gaussian noise that is associated with 
this measurements followed by an optional small letter indi 
cating dimension: d for degrees, h for Hertz and s for 


SECONGCS. 
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The interpretation of the code presented in Figure ee is 
then that both buoys 1 and 2 provide bearing and frequency 
measurements with a period of 100 seconds between consecu- 
tive set of measurements. Buoy 1 provides its first set of 
measurements 100 seconds after the start of the problem; 
buoy er, perhaps because of its limited ranger provides its 
first set of measurements only a ter 500 seconds of the 
start of the problem. The standard deviation of the measure= 
ment noise is 5 degrees for the bearings and 0.04 Hertz for 


the frequencies. 


The filter receives the measurements for processing. It 
iS given some initial conditions and the initial position is 
marked in the figures by the first diamond with the letters 


fee eon top. 


As the target moves and the measurements are simulated, 
the filter generates estimates of the target parameters 
which are represented by smal! diamonds, interconnected by 
dashed lines to represent the estimated path. The diamonds 


and the "+" symbols correspond to the same time points. 


Associated with each estimated point the filter provides 
an estimation error covariance matrix which indicates the 
amount of confidence or, conversely, uncertainty that should 
be given to this estimate. An interesting way to show this 
in a trajectory plot like Figure eed is by plotting the one- 
sigma error ellipses associated with each estimated point. 


These error ellipses have their orientation and axis dimen= 
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SIONS Qiven by Equations (4.27) = (4.29)4 and can be seen in 


Figure 23. 


These figures were obtained by running the’ simulation 
Once. They may represent a typical real-time result. 
Nevertheless they do not give the expected or average 
behavior of the filter to the situation simulated. To obtain 
this average behavior many Monte Carlo runs have to be exe- 
cuted and their individual results processed to generate the 
sample statistics. The sample statistics formulas for 
obtaining the sample means and covariances can be found, for 


example, in page 86 of Ref. [fe]. 


In this work all the statistica) results were obtained 
from 200 MonteCarlo runs. Figure e4 shows the average 
behavior of the filter after these runs. Now the small] dia- 
monds represent the sample mean estimated position and the 
ellipses show the spreading of the results about the mean 
values. Note that at the initial position points since in 


all runs it was the same, the ellipse is the ooint itself. 


To observe the errors in the estimation of the speed and 
heading of the target different plots have to be used. Fig- 
ure 25 shows the sample mean error in the estimation of 
speed as a function of time. Figure e6 shows the sample 


mean error 1n the estimation of the heading of the target. 


Other types of plots or outouts could be generated “from 
the results given by the program simulations, like the sample 


covariance between estimated velocity and bearing, tf ome 
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ever becomes necessary. 


B. SELECTED RESULTS 


As has been seen in this work, the number of variables 
in the passive tracking problem is very largee Some specify 
the true target track and the true position and characteris= 
tics of the sonobuoyss others define the measurement 
scheduling and noise characteristics. The nonlinear filter 
that processes the measurements has its own variables and 
must also assume values for the other variables, which may 


be different from the true ones. 


For any combination of these variables a new problem is 
formed and its simulation may generate a number of interest= 
ing conclusions. The results presented are representative of 
those obtained from the execution of the simulation program 
described in Appendix E. In most of the results presented a 
certain reasonable initial condition for the filter was 
assumeds the position of the buoys was assumed constant and 
knowns the ranges within which reliable bearing, frequency 
and time delay measurements could be obtained were assumed 


equal ° 
1. Nonemaneuvering Target 
ae One Buoy in Action at a Time 


Figure 27 shows the behavior of the filter in a 


Situation where only one buoy 1S 1M contact with the target 


ee 





at any instant. The range and position of the buoys” were 
presarranged to satisfy this requirement. Target speed is 
5.0 m/s? its heading is 315 degrees measured from the posie- 
tive X-axis. The initial assumptions of the filter include 
6.0 m/s for speed and 325 degrees for heading. The measure- 
ments were processed as suggested in (2). Frequency and 
bearing measurements are available every 100 seconds’) and 
orocessed simultaneouslys standard deviations of 5S degrees 


and 0.04 meta 2 were assumed for the notse components. 


The frequency measurement provides, indirectly, 
some range information ands this ways complements the bear- 
INQ measurement. This range information can also be obtained 
by judiciously positioning the next buoys or by having two 


or more buoys im contact with the target at one time. 


In Figure 28 the frequency measurement was elim= 
inated and the simple bearing measurement allowed the estio- 
mator to track the non=maneuverina target almost as wel! as 
before, because of the special position of the buoys. Note 
that the ellipses tend to align with the measurement lines. 
This can be explained from the discussions of Section 1V,B 


and Appendix D. 


There are many ways to improve this tracking. A 
very small improvement can be obtained if less noisy or more 
frequent frequency measurements could be obtained. To get 
more precise frequency itndication one needs better resolu= 


tion from the FFT processors; to obtain better resolution 


128 





fel B. 160 S520 


neo lh as 
i uoy f %*~@O Riv z 


-2.68 Vv a4 > 
eS / 
(km) L e 
eo a 
~4.69 &- : 
° a 
/ & Lr 
Uw, 
/ Yo 
my 
-§.@9 at MS) Y 
voy 3 yo _ uoy 4 
. ; 4, 
-12.¢8 VY’ = 
‘Se 
-12.88 


-4.60 -2.60 0.00 2.06 4.06 6.60 8.40 10.68 12.¢9 
(km) 


® 


Figure 27 : True and estimated trajectories for a 


eee ee eee ee _—_—-—= - 
~- ~- = SS eee eee eee 


non-maneuvering target - Simultaneous 
frequency and bearing measurements. 


129 





¢ B 100 S.@d 


y, 
2col QR 


8.29 
bucy f buoy 2 


-2.66 


(ke) 
-4.80 ¢ [| 
y 


-6.68 


no 
~8.e08 a 3 
uoy 3 a " 
-18.¢e9 


12.60 | &, 


-4,¢93 -2.@0 2.88 2.88 4.€9 6.e9 8.68 16.@8 12.@9 
(ka) 





_..._.__.._.Figure 28_: True and estimated trajectories for a_..._......_. 
non-maneuvering target - Bearing mea- 
Surements only. 


130 





one needs longer records” which, in turns tntroduce more 
errors due to target*sensor geometry variations during 
record time. To reduce the period between consecutive 
independent frequency indications one has to sacrifice reso- 
lution but this would reduce even more the effect of the 
frequency measurement on tracking accuracy. Besidesr an 
emitted frequency of 300 Hertz from a target with a relative 
velocity of 5S knots toward a buoy is shifted by only about 
0.3 Hertz and it can be seen that not much resolution can be 


spared. 


The bearing information is basically continuous 
and it seems reasonable to admit that only economical limie- 
tations exist for obtaining less noisy or, equivalently, 
more frequent bearing measurements. Figure 29 shows the 
effect of processing only bearing measurementsSs but more 
frequently (at eS5-second intervals). In Figure 30 the stan- 
dard deviation of the measurement noise was reduced to 1 
degree. Note that the alignment of the ellipses with the 
measurement lines is more pronounced with less noisy meas- 
urements. Note also that the error in range while the first 
buoy 1S IM contact with the target iS not corrected by 


improving the accuracy of the bearing measurements. 


Figures 31 and 32 Show the mean errors in 
estimating the velocity and bearing of the target. The syme 
bol "+" jis associated with Fig 27, the triangles with Fig 28 


and the diamonds with Fig 29. 
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Since the contribution of the frequency measure- 
ments 1S so smal] compared to that of the bearing indica- 
tions in this case, the partitioning of the measurements has 
very little effect on tracking accuracy. Also the iterative 
techniques suggested in Section IV,B,2 have no marked effect 
on the processing of the frequency measurements for the same 


reasone 


b. Two Buoys in Action at a Time 


The placement of another buoy in action can 
greatly improve hye tracking accuracy, as shown in the next 
figures. The buoys are now arranged so that two buoys, 
instead of one, are in contact with the target at any time. 
The scaling 18 changed so that the estimation errors may be 


better appreciated. 


In Figure 33 the added buoy is a Lofar buoy that 
can only provide frequency measurements. The two frequency 
and one bearing measurements available every 100 seconds 
were processed simultaneously. Note the improvement in accu 
racy when buoys 3 and 4 gain contact with the target in 


replacement of buoys 1 and ed. This also is explained from 


the discussions of Section IV,B. 


In Fiqure 34 a1] buoys contribute with frequency 
and bearing measurements which are still] processed simulo- 
taneously every 100 seconds. The position of the buoySs in a 


line perpendicular to the true track, i8 not avery good 
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position for two Difar buoys. Nevertheless the resulting 


tracking accuracy 1S agood for this non=maneuvering target. 


In Figure 35 the two frequency measurements, 
which require a five=state plant for their processing, were 
replaced by a frequency difference measurement with 0.04 
Hertz of standard deviation in its noise components which 
requires only a foursstate plant, as discussed jin Section 
Il,0,e. The results are basically the same, mainly because 
they are mostly dependent on the bearing measurement and not 
on the frequency information. The use of the frequency ratio 
measurement also gave the same kind of resultS when avery 
low noise measurements, with about 0.0001 units of standard 


deviation, was considered. 


In Figure 36 the frequency measurements were 
eliminated and only Se AphinG measurements were processed. 
Because of the bad position of the buoys, on a line perpen- 
dicular to the true track, the accuracy of the tracking 
along that line is not good. The ellioses show a reasonable 
spreading of the estimates along perpendiculars to the true 
track. Figure 37 shows the same plots, with different scale 


ing, for the first 1000 seconds of the path. 


For different positions of the buoysSs especially 
if one moves the buoys off the perpendicular to the true 
me 
track, better tracking can be obtained as Shown in Figure 


538. The tendency of the ellioses to align with the measure 


ment lines explains the improvement obtained. 
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By reducing the period between consecutive bear-=- 
ing measurements and processing them separately one may 
obtain the accuracy shown itn Figure 39. The ellipses’ are 


like the ones of Figure 38 but with reduced dimensions. 


Good results can also be obtained if one 
includes the time delay measurements. Figure 40 shows the 
result when they are applied to this situation. Notice’ that 
the ellipses are very small reflecting the fact that very 
good information about the position of the target is aiven 
by a low-noise time delay measurement. Also the ellipses 
tend to align with the tangent to the measurement’ lines 
defined in Section IV,B whichs in this case mie teninte delay 


measurements, are hyperbolas. 


ec. Maneuvering Target 


The mext paragraphs discuss results obtained when a 
selected target maneuver is simulated --=- a totai turn of 
180 degrees at 9 degrees per minute with a constant speed of 


SAY) m/Se 


ae One Buoy in Action at a Time ° 


Figure 41 shows the tracking obtained by proo- 
cessing frequency and bearing measurements provided by a 
sinale buoy located close to the center of curvature of the 
target path. The filter takes some time to react to the 


target maneuver ands, when it doeSsr an over=correction is 
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applied. The result is that the estimated path is very dif- 
ferent from the true path of the target. About the same 
result 1s obtained when the buoy is-placed on the other side 
of the track, 1.@.er a small delay to react to the maneuver 


followed by an over-correction. 


In Figure 42 the frequency measurements were 
eliminated and only bearing indications were processed to 
generate an estimated track very close to the one of Figure 
al. The alignment of the ellipses with the measurement 
lines 18 again very clear and may suggest where two buoys 


should be placed in order to improve the tracking. 


be Two Buoys in Action at a Time 


Figure 43 shows what can be obtained by process= 
ing the measurements provided by two buoys. The range of the 
buoys was adjusted so that both are in contact with the tar-~ 
get during all the path shown. Both are Difar bouoys and pro- 
vide bearing and frequency measurements with a period of 100 


secondss, which are processed simultaneously by the filter. 


In Figure 44 the measurements were processed 
separately and a small] improvement in the tracking was 
obtained at the end of the pathr, although during the. 
maneuver the simultaneous processing worked better. It was 
observed in other simulations that the degree of improvement 


obtained by partitioning the bearing and the time delay 


measurements is dependent on the maneuver of the target and 
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the relative position of the buoys. The worst situation is 
when the target maneuvers toward the buoys. This happens 
when the partial estimates are worse than the predicted 
values and occur sometimes during fast maneuvers, or slow 
maneuvers with long intervals between consecutive measure= 


mentse 


In Figure 45 the low=noise time delay measure- 
ment was added and, instead of reducing the errors, the 
tracking became worse. The direction of the ellipses is 
clearly along the tangent to the hyperbolas of position, 
indicating that the filter is iiagat cain about the range of 
the target to the buoys during all the maneuver. This 
characteristic of the measurement lines in this cases, and 
the fact that the target is modelled by the filter as ee: 
lowing basically a straight path, explain the form of the 


estimated track. As soon as the target reassumes a constant 


path the filter starts to recover. 


In Figures 46 and 47 the iterative techniques 
described in Section IV,B,2,b were applied to the frequency 
and time delay measurements. In Figure 46 only one itera- 
tion was executed; in Figure 4/7 enters ted formulas were 


applied three times. The imorovement is clearly seen by come- 


paring with Figure 45S. 


As with the partitioning of the measurements, it 
was obvoserved from other similar simulations that the applio- 


cation of the iterative techniques does not always provide 
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better tracking. It depends on the type and direction of the 
maneuvers, the position of the buoys and, consequently, on 


the position of the measurement lines. 


In the case of the non=-maneuvering target the 
use of two buoys instead of only one was almost 5 nance suf- 
ficient to provide qood tracking. If for the maneuvering 
target one uses three buoys judiciously positioned, instead 
of only one or twos one can also obtain good tracking as 


clearly shown in Figure 48, 


C. SUMMARY 


Below are listed the most important facts observed from 


the simulations: 


“== The tendency of the error ellipses to align with the 
measurement lines was observed in Section IV,8. The simula=- 
tions show that the sample covariance ellipses also follow 


this trend. 


“--=- The tendency of the error ellipses to align with the 
measurement lines can be of oagreat help in the practical 
solution of filtering situations. If. for example, the error 
ellipse is very thins i-«@e,r have a very high ratio a/b, then 
the best measurement to process is the one whose measurement 
line is approximately perpendicular to the major axis of the 
ellipse. In the trackina oroblem this can be obtained by a 


bearing measurement from a buoy placed along the perpendicu= 
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lar to the major axis of the ellipse, or by a time delay 
measurement from two buoys placed along the direction of the 


major axis of the ellipse. 


~-- The frequency measurement have a very small effect 


on the tracking accuracy. 


“~-- The use of frequency difference or frequency ratio 
measurements improve computing efficiency (by eliminating 
the need for having the rest frequency as a state in the 


filter model) without noticeable effect on accuracy. 


w-- The partitioning of the measurements improves come 
puting efficiency, may provide better tracking accuracy and 
allows the processing of measurements as they occur, and 


thus is of great practical importance. 


--- For low noise measurements the iterative techniques, 
although requiring an increase in computinas, are capable of 
providing better tracking for maneuvering targets. 

--- Tracking with only one buoy is acceptable only § for 

\ 
non=maneuvering targets. 

--=- The relative position of the buoys is one of the 
most important factors which influences the quality of the 


tracking. 
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VI. CONCLUSIONS 


A. SUMMARY OF RESULTS 


The optimal estimation of characteristics and/or parame= 
ters of a certain class of nonlineary, dynamics stochastic 
systems has been studied in a probabilistic environment 
using Bayes formulation concepts. Approximate solutions and 
filtering algorithms were generated ands, amonq them, the 
known Extended Kalman Filter equations’ and higher order 


filtering equations have been obtained through this method. 


The problem of tracking submarine targets using special 
passive sonobuoys was modelled and with this model extensive 
simulations were executed to allow the study of the problem 


in detail. 


Most of the results indicate that the frequency measuree- 
ments have minimal effect on the filtering process. The 
small contribution to range information that thev do pro- 
vide, when associated with bearing measurements, can nor 
mally be obtained otherwise by judicious placement of subseo- 


quent buoys. 


The utilization of other types of measurements Such as 
the frequency difference, the frequency ratio and the time 
delay measurements, proves to be a great help in improving 


computing efficiency by eliminating the rest frequency as a 
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necessary state, and thus reducing the dimensionality by 
onee AlSoOr especially with the time delay measurement, a 


great Improvement in tracking accuracy is possible. 


The concept of partitioning the measurements and pro- 
cessing them separately, even if they occur simultaneously, 
1S Shown to bring advantages in computing efficiency and 
alsor, for nonlinear measurements, In tracking accuracy. 
This concept is of great practical mmoortance, forr In 
Situations where the measurements are sporadically and none- 
periodically received it 1S most important to be able to 


process the measurements as they naturally occur. 


The graphical interpretation of the action of Kalman 
filters, developed in this works provides insight into the 
importance of each variable of the problem in the filtering 
process. The direction and magnitude of the correction which 
1S applied to the predicted values to generate the new estico 
mate values can now be anticipated as a function of the 


measurement. 


Nonlinearity errors have been graphically presented and 
iterative techniquess including the known Iterated Extended 
Kalman Filter equations, have been suggested to counteract 
their disruptive effect. The application of these tech-= 
niques to the tracking problem shows that improvement in 


tracking accuracy 1S possible. 


The granohical interpretation also indicates the very 


practical conclusion that the error ellipses tend to align 
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with the measurement lines, as defined in Chapter IV This 
provides guidance for optimal positioning of the buoys and 
the types of measurements to process. This observation 18s 


reinforced by the simulation results of Chapter V,. 


B. SUGGESTIONS FOR FUTURE RESEARCH 


The concept of partitioning the measurements allows one 
to. process the measurements Separately and opens some ques= 
tions for future research, such aS-e in which order should 
nonlinear measurements be processed to obtain maximum effi- 


ciency? 


The graphical interpretation of the filtering process 
allows anticipation of the direction and magnitude of the 
corrections which are then applied to predictions to gen=- 
erate new estimates. This situation suggests future research 
in the determination of the optimum characteristics of the 
measurement functions which in turn can determine the 


optimum positions and characteristics of sensors. 


The expansion of the model to include accelerations 
should be considered if maneuvering targets have to be 
tracked efficiently and the increase in computing power can 
be obtained. Consideration of the depth of the target and 
the inclusion of uncertainty about the vosition of the buoys 


are additional ways to extend this study. 
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The model and simulation implemented herein, without 
benefit of classified information, is very flexible, but 
many idealizing assumptions have been made. Ihus, a natural 
extension of this work is to apoly the algorithms and con- 


cepts developed to a real problem using actual sensor data. 
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APPENDIX A 3; 


The relation 
pee te aN 


1S given where x and 
and qr respectively, 


tor and G is a (nxq) 


PROBABILITY RELATIONS 


Cares) 


v are random vectors of dimensions on 


=— 


ais a n-dimensional deterministic vec- 


matrix of deterministic coefficients. 


The joint density p(y) is agiven and the joint density 


p(x) is wanted. 


Case (i) - If 
given by (15), 

v= uate - aj = 
and 

ex). = joer| a 


Case (i131) - If 


Sori “nti” ~~ 


and create 


f 


= 


n+2 — 


Nn = q and Gi existS, the solution 1S 
f(x) (Age) 
| | - o (Ff(x)) (A.3) 

a 


nm < q one could add dummy variables 


ren) I eee? *q = ‘q 
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where 


Now, 
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lw 


|x 
' 
‘ 
| 
t 

lm 
| 


jo 


i f Ca exists, 


i 6 {xt = a} = t(x"™) 


and 


apoli 


Let 


From Case (i), the solution is 


') = aoe 
aux) = joe 4] | 2 (CEG) ) 


esi ei) Sit n > Q the following 


eSe 
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ee 8 trey 
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[x*' - a‘) = f(x") (A.7) 


From Case (1) one then has 
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a Wp at Em 
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Bee 7 eG? = 





qt1 
© ORIG? 5; TOIL en = STN ela sta 
but 
Seti rey, Seesiehe  Shl yO 29 Slangin @ 
- Sara Sail ea 
or 


165 





-_ ? ~- 


From (A.9) comes 





i) - - 
oO a UG 2 2 OG OU 7 - SE, VeHine 
and 
p(x Xx GFeeos X ) = 6 (x = } 0 (x ax peeoeg X ) 
2 qt+1 qt+1 qt+1 me Gell sz q 
Doing the same thing for x Peoria 4: , one finally 
qt+2 n 
gets 
n 
p (x) = 6(x. © 1.) . set| | | oe Cx. 2 
al ale ax yV—-c- 
ier (A.10) 
moe the special case of gq = 17 the vectors x, (6) and 


g' become scalars and the sbove equations lead to 


‘ = .¢ = = / 
f(x") (x)? Cx, ae a 








1/aq é (C(x 2a )laq ) 
1 ny 1 1 ar 


i=2 (Moll) 


Another similar formula can be obtained for this spree 
cial qlC« = situation. The matrix Gin (A.1) 318 Now a ne 


dimensional vector and, if one breaks (A.1) into its on 
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individual scalar equations one can obtain 


 ] = = = 
p(x iv) = 6 (x, a5 gv) 


Also note that 


and 


Bx o9X 1V¥) =wol(xo,x.rVv) p(x, tv we o(x, ivi.elx iv) 
L Al 1 Jj J x J 


And from these relations,s 
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APPENDIX B =: COMPONENTS OF COVARIANCE MATRIX 


For a single measurement of the form z = H.x + ve the 


—— 


observation matrix becomes the vector 


where n 1S the dimension of the system. 


The Kalman gain is given by 


T 
G = P'H (HP'H! + F7} (B.1) 
where P' is used to represent the prediction error covari-= 
ance matrix, which 18 symmetric, and r is the standard devieg- 


ation of the meaSurement error. 


Since H is 1M thiS case a vector, the product HP tHE 


will be a scalar. 
Let c = CHPT'HD + p*7} (pace 


The vector P'H can be shown to be 


Stee 
jae es 
" | 

» a} 


2%»; : 
P'Ho =| jal “J 4 (B.3) 


9 Saget oe 
1 8 J 


imps 


j 


and thus the constant c 1S given by 
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n n 
) ee Pa yamene in) ec! + rj 
i=l j=l ij i 3 


Cas (B.4) 


From (B.1) one now sees that the Kalman gain 1s a vec= 


tor whose components are the components of P'H! multiplied 
by the constant Cr 1.e. the components of the gain vector 
are 
n oe 0) a) 
a = ja P| h. / (2, oil Pi ye itis vip) (62>) 


The estimation error covariance matrix is computed from 


the equation 


mee tl = G6. Hi eie= Cee (26) 
where 
AP* = = GHP' = = cP'H HP' (B.7) 
Let's define 
a, = e See =p! ah. + 6° heart cacao A (B.8) 
if j=l a Ee ae dab oi. 2 eZ in n 
Then, from (B.3), 
HP' = Touro) e = fa 3a eden ama (B29) 
. 1 2 n 
and from (B.4), 
n 
2 bey ce veer po) (B.10) 


From (B.7) the tincrement 


matrix 


AP’ to the 


is then 
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APPENDIX C s CHARACTERISTICS OF MEASUREMENT FUNCTIONS 


Ln Figure C=! some curves of the form f(x) = constant 


are shown, where f(x) is assumed given. The vector x* and a 
b 


constant value d are also given. The vector x is defined 
by the equation 
qd h 
Ee a ae (C-1) 
Ie} dT 
or 
Ros x? + he td7 hua (een 


where h is the gradient of f(x) with respect to x taken at 


point i xe and h/thi is a unit vector along h. 


. 8 ye 
y= ees + hoe: bak = Coan Le 
ee Tes Wes a 
al OX) ~ y dX, =~ 


lo determine the location of x numerical data 1s 
required. However, it is possible to determine qualitatively 
whether xb ends on the curve f(x*) + dr or on one side or 
the other of this curve. This is done in the following 


paragraphs for three special functions. 


b 


Case (i) For linear f€x)d, x always ends on the curve 


f(x) = (3) + dy as shown below 





x) = f (x°) asl 


f(x) = f(x )~_ 


E(x) = f(x) ‘ 


b |a| / {hy “ 





Figure C - 1 : Problem geometry. 
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Let f(x) = px. + qx 


V 2 


Ihen h, =p ’ ho = q and | n |? = 0° ¢ of 


From Equation (C2), 


x? = x* + (pe, + ge, ).a/(p* vs q*) 
So 
a8 = nt + ed/(p2 + q*) 
i = a + ad/(p2 t q*) 
Now, 
f(x?) = px + axe = px, + p* a/(p2 + qq?) + 


+ ax? + at d/(p* + qe ) = 


TT 

se] 

bad 
hr 


+ ax + ds 1x) tC 


b 


Case (i1) = For f(x) defined as below, x always ends 


on top of f(x") + d, as in the Case (j). 


(x) = Ux = p)® + (x. = 34) 1/2 
Now, 

es 6X4) 7 PIT m ae OS Se ony 
where 

m= exe - p) + (x om q)2j1/2= mex) 
and 

SA 7 ei 

oh a + nS l 





From (Orc )e 


b . a . a aoe a. 

x Sx + hed = x ¢ CCx, pie, + (x, aje J ~d/m 
So 

eax? fohiice=) ounce enya 

‘l 1 1 1 il 

baa) =a - a a 

Ky = Xy + h,d = x, + Cx, qa)d/m 
NOws 


f (x is ((x? + hid = p) + OF t pe = q)2 V2 
‘ a 2 a be 

rm en - + h + 
E(x) p) GeGx p) d 


t (x2 = g)? + 2h d(x? - foe 
(x, q) oe q) ; } 


Using the values of hs and hoe 


#(x?) = UCx9 - 5) + (5 - ay + g? + 


+ 2d UC x3 - p) /m + (x3 ~ q) /mj)/22 


f(x? + 2d.fCe2) + a 2yt/2e 


1/2. 


ioe 0 


FOL? + d 


Case (iii) = For f(x) as defined below and shown In 
Figure Cec, x? always ends before the curve f(x) + ds for 


[a| <1. 


f(x) = arctanl(x, - q)/(x, = 9)] 
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Figure C - 2 
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lhe case la| > t/e is not of interest since the lines h 


and Pox) = axa) + d are diverging. For la| <n /ep 


a = © (x5 = q)/n h, = (xt - p)/n 


where 


n= (xt - p)? - (x5 - q) 


and 
2 Z 2 
= + = J 
|r hi ho l/n 
rrom Equation (Cee), 
x> = x + de (= (x5 - qje, + (x) = De 
SOvr 
| a: ae a. 
x, = xT (xe a)d 
ne = x5 + (xt - p)d 
Nows 
(x2) = arctan{(x? = a)/(xP = pil = 
= arctan ([(x4? = q) + (x? |= p)d)/ 
ae 1 
(C(x? = po) = (@ = q)d) 
1 iayZ 
or 
t (x? ) = arctanl(s + d)/(1 =- sd)] 
where 


Sue (x) - q)/ Cx - p) = tanlf(x?)) 





SO, 


cant (x’)) = (s + d)/(1 = sd) = 


[tanif(x )) # d} 7 (1 = dotamlt (x ))10- 


tan(f(x ) + arctan(d)]) 


t (x?) = f(x" ) + arctan(d) 


Since Lal < m/e, then arctan(d) < d= and xb will not 


reach the line f (x %) + d. 





APPENDIX D s ANALYTICAL EVALUATION OF KALMAN FILTERS 


After the design of a filter an analysis phase is 
necessary to verify its behavior in various representative 


Situations. 


A "filtering Situation” is considered here to be come 
pletely ee wan fice by: 
—aoetnemerdesperametersmand initial conditions of the 
system which Generate a uni quem track x (0) eG), sa. 7 0 tm 
“-- the true parameters of the sensors and their meas= 
urement schedules. 


“<= the parameters and initial conditions assumed by 


the filter for the system and for the sensors. 


Ihe approach normally used to determine filter behavior 
is to simulate the desired situation and execute hundreds or 
thousands of Monte Carlo runs and compute sample statistics. 


The most useful results of thisS process are: 


a (k) - the sample mean of the estimated state vector 


at time Se obtained from m Monte Carlo runs. 


e tk) = a (k) - x(k) = the sample mean of the 


estimated error vector at time tie 


b (k) * the sample second central moment of the state 


(or error) estimates about A Oe at the 
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the objective of this section is to study the possibil-=- 
ity of obtaining values of alk), e(k) and b(k), which are 
approximately the values of a (k), em(k) and bu (k) when m 15 
very large, without the use of the time consuming Monte 


Carlo operations. 


1 -"- Theoretical Solution 


Unce a “filtering Situation", as stated above, is 
defined, the operation of a Kalman Gilites is described by 
the recursive equations below where only the dependence of 


each term on the previous estimates is stressed. 


X(ktlik) = f(x k)) (D.1) 


PCktLik) = g(x(k)).PCkK) go T(K(k)) + 


t 9((k)) O(k) 2g) (K(k) (D.2) 


G(x(ktlpk),PCktl$k)) = PCK+LIK) oH Cx(ktlik)). 
(HC2Ck+LEK)) oPCKALIK) HECK CKELIK)) + ROCKIT 


(D.3) 


K(ktL) = xCktlik) + GORCKHLIEKI,PCktlik)).thlx(ktl)) + 


+ vektl) = ACx(ktlik))) (D.4) 
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PCkt1) = C1 = GOxCktlik) PCkeltk)) HOXCkK tL ik) )) PCktlik) 


(D.5) 


In this case these equations describe a deterministic 
process were it not for the random measurement noise vl(ktl), 
the only random input to the "filter dynamics" since a 
unique track iS considered. The initial condition of the 


filter, x(0) and P(0) are also deterministically given. 


tquations (0.4) and (D.5) can be rewritten as 
XCktl) = oCxCktlLikdePCkttik)) + vy CKCktLI kK) /PCk+lik)) ov(kel) 


P(k+l) = FCx(k+lik)»P(k+lik)) 


where $¢ +» ¥ and TI are generally nonlinear matrix functions 


of the variables indicated. 
Considering (D.1) and (D.2) one can then write 
xCktl) = fF Cele ePCkK)) + O*Cxlk) PC kK) oviktl)  (D.6) 


P(k+1) = h*¥(x(k),P(k)) (O27) 


‘where hy g*, h* 


are also generally nonlinear matrix func 
tions. From (0.6) and (D./) one can clearly see now that, 
if v is a discrete white Gaussian noise process, the joint 


process {x,P} is Markov. 


lf one considers a new vector y formed with the n ele- 


ments of x and the n.e(nt+1l)/e distinct elements of the (nxn) 
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Symmetric matrix Py, y is an [ne (nt+3)/2]l -dimensional process. 


Equations (D.6) and (D.7) can then:be combined into 


y(k+1) = ACy(k)) + Bly(k))ev(ktl) (D.8) 


where A IS a vector function and B a R to Rm matrix funce 


tions, m = n.(nt35)/e. 


The compJete behavior of the filter can be described by 


the probability law of yr which is obtained from 


ply(k+tl)) = ol(y(ktidiy(k)) ply (k)) dy (k) (D.9) 


and ply(ktley(k)) is obtained from plv(kt1)) and Equation 
(D.8) in the manner shown in Appendix A. The initial value 
y (0) is deterministic since it contains the initial condi- 


tion of the filter which is given. 


c --- The Linear Case 


For the special case of linear dynamics and observa- 
tions, the solution can be obtained in a simpler way. In 
this case the functions +7 9 and H are not functions of the 


state estimates and so neither is G. 
Equations (D.1) = (0.5) can be grouped now into 


KCKtL) = OCKCKHLEKI PC kt ed) # ¥ (PC KLE K)) ev (kel) 


PCk+1) = f (PCktlik)) 
and finallys in a recursive way, 
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XCkt1) = #°CKCk)/P(K)) + g (PC K)) ov Ck) (Oran 


P(kt1) = h(P(k)) (D.11) 


From (0.11) it is clear that the estimation error 
covariance matrix 1S now a deterministic process and can be 


precomputed, together with the gains. 
Equation (D.10) then becomes 
KCktl) = F(R) + gh (kel) ev ( kt) 
From (D.1) = (0.5) the values of ¢ and g* can be found 
to be, for this linear case, 
KCktL) = CL + Glktt) HC kt1)) 66 (kK) ox Ck) + 
+ GCktl)HCktl) .KCktl) + GOktl).vktl) 


or 


x(k+1) SC kt1).x(k) + F(k+1) + GCktl).v(k+l) COe le} 
where 


S(k) = C1 + G(k).H(K)I.6(k) 


FCk) = G(k) HC kK) ox (kK) 


ihe sought after behavior of the filter can be 


described by the moments below, which agree with [12] 
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a(k+}) = EtxGk+1))) Ss SiG@is 1) 5c (eis eee 


+ GCk+t1).E lv(kt1)] (0.13) 


e(ktl) = a(ktl) = x(ktl) = S(k#l).EEx(k)) + 
+ (1 = GCk+tl)sHCktl di .xCktl) + 
+ GCktl).Elv(k+i)) = 
= S(k+l).e€k) + DCktt).x(ktl) - SCktid.x(k) ¢ 


+G(kti) Elvktt)) (D214) 


D(ktl) = S(ktl).b(k).S7(kt1) + 


+ GCk+1) Elv(k+1) evi (kt1)IGE(k+1) = (0.15) 


3 *-- General Case { 


For the general nonlinear problem one returns to Equa- 


tion (0.9) which cannot normally be solved in closed form. 


Numerical techniques and approximations would be necese- 
Sary that may use more computing power and present worse 
results than the Monte-Carlo process that we are trying to 


avoid. 


Une way to simplify the problem is to assume that y(k) 
has an approximate Gaussian distribution, even with all the 
nonlinearities shown. Linearizing Equation (0.8) around the 


mean value of y(k) would give 
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y(kt1) = ACy(k)) + - ACy(k)).ly(k) = y(k)) + 


+ Bly(k)) ev (ktl) (Okie 


Taking y(0) as the deterministically known initial 
value yO), and applying the expectation and variance opera= 
tors to Equation (D.16), one gets the recursive equations 


for the moments of yk): 


y(k+1) = ACy(k)) + Bly(k)) Elv(k+1)} (D.17) 


Varty(kti)) = C2 AG(K))) Varly(k)}.f © ACy(k))I + 
Seles = ie = 


+ BCy(k)) Variv(ktl)) BT (y(k)) (D.18) 


and 


ha 


jhe moments we are interested in, 


| | 
lol 


v are 


directly obtained from subvectors and submatrices of the 


above moments. 


The primary difficulty with this aporoachr however, 15 
in obtaining the functions A and 8B. This can be seen for the 


very simple case belowr aS an example. 


Suppose a scalar system with a single observation 1s 


described by the equations 
eC k1) = x2(k) + wk) 


z(k) = x(k) + v(k) 
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Equations (D.1) = (0g5).46 1.4 


RC ktLik) = x 2Ck) ; e(k) = 2x(k) 
ay) 2 

p(ktlik) = 4x Chui k ice O null) 
WEE (a) opal ss + a7 (k) 

G(ktl) = 


—_— 2 ~ 


Pee Napa) ee Gea) Ge a2 


and, after the appropriate steps, 


; . ux Ck) ep lk) + a7(k) 
ck +t) he 
, ax“(k)eo(k) + Q°(k) + r“C(k) 


(ux? (k) p(k) + Q*Ck)) er? (k) 
p( k+l) ae eee 
_ x ux-(k)ep(k) + a2(k) + 92k) 


4x7(k).p(k) + Q7(k) 


WaCun ean) @ ele moa 
« ev( k+l) 


The mean value given by Equation (D.!7) can be obtained 
in a simpler ways however. It is the result obtained by 
running the filter in the simulated environment but making 


the measurement noise assume its mean value Elv(k)]+, nore 


mally Zeroe 


a 


If one uses more terms in the expansion of Equation 
(0.8), better results willwbe obtained at the expense of 


increased computing effort and time. 
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For this general case it 1s suggested that these equa- 
tions be applied to a simole scalar or twordimensional none 
linear problem and that the computing power required to find 
the sought after moments be compared with that required to 
run a Monte Carlo process yielding the same level of accu- 


racy in the results. 
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APPENDIX E : COMPUTER PROGRAM 


Ihe basic structure of the computer program written for 
the simulation of the tracking problem, and the evaluation 
of models and filtering algorithms, is shown in Figure Eel. 
The solid lines represent the normal flow of control within 
the program, the dashed lines show the extra pathes’ that 
become available whenever an interruption is requested by 


the user. 


After a brief introduction to the program a table 
called MENU is presented to the user, as Shown in Figure Ee 
ele if the user chooses actions number 1 or 5 the result is 


clear. 


Choice number 3 provides access to all the problem 
variables. Since the number of variables that characterize 
each simulation 1s very large, a simple way to deal with 
these variables had to be devised. As it 1s implemented, 
the program always "remember" the values given to the varie 
ables at the last time the program was usedr so that only 
the variables to have their values modified have to be 
addressed. The change of value of a variable is simply made 
by selecting the page where it appears, typing the letter 
associated with the-variable and the new desired value. The 
program immediately responds by presenting back the entered 


value. A more detailed diagram of the actions resulting 


~ 


187 





Introduction 


Tutorial MEN U . STG Ve 


Information 





Chanae __iChanae Value 
progmam | 0 mor) fone m | 
Charact. | Variables ! 





CPC 


_ wg 


Simulation | 





Printing and 
= tee & J (alle, 
Routines 







Praqure Ewll*: Bastcestetict teem ot OrOogran. 





MENU 


Press the numvoer corresoonding to the action 


desitred and c/r. 


t = PRESENT TUTOREAL INFORMATION 


2 - MONTIFY PROGRAM FLAGS 


3 © FORMULATE OR MODTFY THE PROBLEM 


go = SI4kt THE SOU OIG Cl THe PRG bbe 


S - FENN THE PROGRAM ANN EXTT 


FPyqure E=c < [me me Nletaal >. 


So 





from choice number 3 is presented in Figure E=3. A typical 


page of variables is shown in Figure Ed. 


Choice number 4 simulates the tracking problem accord- 
ing to the values defined in the previous step. A block 
diagram of the actions involved in the simulation is 


presented in Figures E=5S, E-6 and Ewr7. 


Initially all the important problem variables’ are 
printed with their current values in a form as shown in Fige- 
ures E=6, E=9 and E-10. During the first run the scheduling 
of all the events involved in the simulation is also printed 


for future references as shown in Figure Ew ll. 


With choice number ec from the MENU anew table is 
presented as shown in Figure Eel2d. This table in also 
presented to the user whenever he requests an interruption, 


at any time or point within the program. 


The extra choices that now become available are almost 
self-explanatory but it should be added that with choice 
number 7 the simulation is ended and the partial statistics 
are computed and written in the appropriate output files; 
choice number 6 allows the modification of any variable in 
the middie of the simulation, in the same way as explained 
beforer choice number 0 transfer control back to MENU or to 


the point of interruption. 


A listing of the program, in C language, is available 


at the Electrical Engineering Department, Naval Postgraduate 
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from 
MENU 


from 
EEC 


Present 
First Paae 


of Variahles 




















Choose Action 








Select Print Paae Print Paae 
Variable Tf Jer Was Tf [Jt Was 
To Be Modified Modified 
Modified 

Chanae | Present 


Next Page 
of Variables 


Value of 
Selected 
Variable 


E10) 
ME ANIL 


—at —— ome ——- oe a a ere Ses = — 


CPE 


Fiaure F-3 : Charginn the value of problem variables 


Boe 





To change the value of any variables type the indicated 
letter, at least one snace, the desired value and c/r. 
To see the other variables of the problem, tvyoe 1 clr. 


To go back to menur tyoe 0 c/r. 


( 1 meter = 1.994 yards; 1 meter/sec = 1.94 knots ) 


TARGFT INITIAL PARAMETFRS 


initial x = nosition 3 0.9 km a 
y= postition = 0.0 km b 
speed : Le see Cc 
headina Seo. Uae d 
frequency ee Oa HORe e 


Standard deviation of forcina functions 


speed/sec : 09.9 =m/sec/sec f 
headina/sec ; 9.9 degq/sec gq 
freaquency/sec °¢ Oer0 hertz2/sec h 


Figure E-4 : Tvpical oagqe of variables. 
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from MFNU 









Print Value of 
All Variables 









Initialize Variables 
and Random Number Gen. 
for a New Run 









Print Table 
of Events 











k 


Is taroet to be updated y i ttodate 


- PeocaVe Gers 


pee 
N 
eo 
‘s, 7 
Nee 
Are There Anv y a 
Other Fvents Se bes, 
. ae j els Pe 
Occurrina 2) See he 
This Time? 7 ee 
N 


Figure E-5 : Simulation T 
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time =+ | 


Is time > duration 
Ot eed RU howe ee gf a 


run =+ ]j 


Is run > number pe {5 ) 
_— 
Oy pl inics Gsm 


Compute | 
Sieaot 1 Sitres 


! 
Send Results! 
To OQutout 
Files 


ge, 
MENU 


Figure E=6 = Simulation ti 
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Select Action 
From Table 
Of Events 










Chanae Taraet 
Parameters to 
Simulate 

Maneuver 










Store Values 
of Parameters 
and Fstimates 
for Statistical 
Purposes 


Ts tarmet to be 


updated hefore —— eu ear me OOo ae 
measurements? — “ Target 


Parameters 


Simulate 
Selected 
Measurement 


Netermine 
Predicted 
State 


oo.” . nt 
; Linearize 
About 

| Predicted 
, State 


















[Process 


Measurement 
To Oktatin 
jFstinate 


Is the requested 
numoer of tterations 
satisfied? 2 


It inearize | 


‘About Partiall 


jEstimate | 


Fiaure Eo? : Simulation TIl 
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TARGET INL TIAL PARAMETERS 


initial x = Mmosition§: O20: km 
Y= POS ito m ee: 9.0 km 
speed : 720) 5 m/sec 
heading * 40520) dea 
frequency s SOOLR heetwz 


Standard deviation of forcina functtons 


soeed/sec 
headina/ssec 
frequencv/sec 


0.0 m/sec/sec 
0.0 dea/sec 
O20 Tehert2/ sec 


e@ se ee 


GENERAL NIRECTIONS 


numoer of 
numoer of 


number of 
number of 
duration o 


number of key ooints 


seed 


PIFAR bhucys. 
LOFAR bucys 


3 buoys 
0 buOVS 


maneuvers : ee oman 
runs : en runs 
f runs : C000 sec 


20 points 


: 13339 


average sound velocity: 1500.0 m/sec 
range of the buoys : 8.0 km 


Plier tT Tae PARAMETERS 


inttial x = nositian 3: 27.25 km 
y 7 nosifian 3 20 k mm 
coeed : 6.09 m/sec 
keaaina SS oye 
freiatienc v SU. hertz 


Standserd sevireation of forcine funNcrTiaons 
speed/sec © 0.905 m/sec/sec 
heaaina/sec . (vot degq7vy sec 
frequency/sec : 0.910 hertz/sec 


Fiaure E& 


Listina of oroblem variables 


Eg'6 


T 


V7 aAndt ya 


noan T ®w 


2m 18) 


> 2 7 oD Do 


_ 





fr 


oo coo o 


f 
0 
0 
OF 
0 
1 


INTTTAL COVAPIANCE MATRIX 
x pos y POs sceed eaaa 
x pas 12,0573 i Sano O10 SG Ore Oa 
y oos * Ure 0.0 bh O20 
speed * k 2.0 m 0.0 n 
heada x * 100.0 5s 
freq * * * 
(km, m/sece dearees and hertz crossmultinonlied) 
MEASUREMENTS SCHEDULE | 
huoy tyoe start period mean stdev 
1 B 100 100 0.90 SU Mabeadet 
j F 190 100 O20 90.94 ahijsk) 
Vag a 500 100 0.0 5.0 mnonoar 
Po F 500 100 Be 09.04 stuvwx 
(start and period in seconds) 

(mean and std dev in degrees, hertz or seconds) 
(period must be zero if measurement is not used) 
MENSURE MEIN omoe he pUien we I] 

buoy tyre start neriod mean stdev 

le T S00 1Q0 0.0 9.01 abcdef 

13 T 590 1gOn 0.0 OL0) eh 

ec 3 i Sire 1OO0 0.0 0.901 mnopar 

3 B S90 160 0.0 50 stuvwx 
(start and period in seconds) 

(mean and std dev in degrees, hertz or seconds) 
(period must oe zero if measurement is nat used) 


Fiaqure E79 3: Listina 
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of oroblem variables II 


< wh? OO ew 





TARGET MANFUVERS 


Ist maneuver ! time = 600 sec a 
vdot = 0.0 m/sec/min b 
hdot = 29.0 deg/min G 
fdot = 9.90 hertz/min oA 
TARGET MANFUVERS 
end maneuver : time = 13209 sec a 
vdot = 0.9 m/sec/min b 
ndot = 9.0 deg/min c 
faqot = 0.9 hertz/min a 
RUDNYS PARAMF TERS 
First huoy : tvoe = ODIFAR a 
x = position = 0.0 ker lh 
y = position = -3.0 km ‘e 
Vt Deak: 
bearing error mean= 9.90 dea d 
std dev = a0 deg e 
BUOYS PARAMFTFRS 
ena buoy : tyvoe = OJFAR a 
x = position = Ler), km bh 
y = pasition = -3.0 km Cc 
If ODTFAR?: 
bearing error mean= 9.9 aea d 
Std gaev = S20 deg e 


Fiaqure E-10 : Listina of poroblem varitanles ITI 
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700 


TUNG 


700. 


R00 
ROO 
ROO 
R00 
R00 
R00 
900 
9090 
900 
900 
SNe 


900 


SEQUENCE OF EVENTS FROM 60! [CU 900 SGEC ants 


SEC 
sec 
sec 
sec 
sec 
sec 
sec 
sec 
sec 
sec 
sec 
sec 
sec 
sec 
sec 
sec 


sec 


bearina measurement hy buoy number 


I 


frequency measurement by huoy number } 


hearina measurement by buov number e 


frequency measurement by duoy number e 


time delay measurement bv buoys 


Monte Carlo roint number 7 


beartna measurerent hy buov number 


1 


and 2 


1 


frequency measurement hy ouoy number I! 


bearina measurement by buoy number e 


frequency measurement by buov number e 


time delay measurement by bhucys 


Monte Carlo coint number 7 


hearina measurement hy buov number 


} 


andme 


1 


frequency measurement by buov nunber | 


hearina reasurement hy buoy number 2 


frequency measurement hy bnuov number ec 


time delay neasurement by hunys 


Monte Carlo coint number y 


Fiaure E-11 ¢: Tahle of events. 


ed, 


1 


ania 





PRONGPAM CHARACTERISTICS 
Press the number corresnondina to the option or 


action desired and c/r. 


To continue or return oress 0 c/r. 


1 - update tarcet every second 


e - undate tarcet only nefore measurements 


3 = start orinting on line orinter 


4G = ston printina 


5 - put results of next run on outout files 


6 - modify parameters of the problem 


7 - terminate the oroblem 


Figure E=l2 *eiat cerues omega 6. 
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School, to any interested reader. 
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